{'fraction_non_alphanumeric': 0.03395843892549417, 'fraction_numerical': 0.019656669384517066, 'mean_word_length': 5.274889380530974, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mapping germ-layer specification preventing genes in hPSCs via genome-scale CRISPR screening

Summary Understanding the biological processes that determine the entry of three germ layers of human pluripotent stem cells (hPSCs) is a central question in developmental and stem cell biology. Here, we genetically engineered hPSCs with the germ layer reporter and inducible CRISPR/Cas9 knockout system, and a genome-scale screening was performed to define pathways restricting germ layer specification. Genes clustered in the key biological processes, including embryonic development, mRNA processing, metabolism, and epigenetic regulation, were centered in the governance of pluripotency and lineage development. Other than typical pluripotent transcription factors and signaling molecules, loss of function of mesendodermal specifiers resulted in advanced neuroectodermal differentiation, given their inter-germ layer antagonizing effect. Regarding the epigenetic superfamily, microRNAs enriched in hPSCs showed clear germ layer-targeting specificity. The cholesterol synthesis pathway maintained hPSCs via retardation of neuroectoderm specification. Thus, in this study, we identified a full landscape of genetic wiring and biological processes that control hPSC self-renewal and trilineage specification.


Introduction
Human pluripotent stem cells (hPSCs), including human embryonic stem cells (hESCs) and human induced pluripotent stem cells (hiPSCs), can self-renew and differentiate into all three germ layers (Thomson et al., 1998;Takahashi et al., 2007). Thus, hPSCs can be considered as highly advantageous cells for studying human development, modeling inherited disorders, and replenishing degenerated tissues and cells (Soldner and Jaenisch, 2018). Sophisticated protocols have been developed to target hPSCs into the ectoderm (Ec), mesoderm (Me), or endoderm (En) lineages, primarily via cocktailing developmentally related signaling cues. Activin/Nodal and fibroblast growth factor (FGF) signaling are two dominant pathways that maintain hPSCs in an undifferentiated state, and interference with the signaling molecules of these two pathways will lead to spontaneous differentiation (Vallier et al., 2005;Xu et al., 2005). Once leaving the culture conditions favors pluripotency maintenance, hPSCs adopt a neuroectoderm (nEc) fate in the absence of extra signaling activators, and dual inhibition of bone morphogenetic protein (BMP) and transforming growth factor b (TGF-b) signaling further strengthens this nEc fate (Chambers et al., 2009;Xu et al., 2005). On the other hand, activation of the BMP/Wnt/Activin/Nodal/FGF signaling pathways is major driving forces for generating primitive streak-like cells and Me or En lineage (D'Amour et al., 2005;Tam and Loebel, 2007;Zhang et al., 2008;Martyn et al., 2018;Wang and Chen, 2016).
Extracellular signaling pathways govern hPSC fates largely through their convergence on regulating fate-determining transcription factors (TFs) (Xu et al., 2008). OCT4 (POU5F1), SOX2, and NANOG are well-characterized core pluripotent genes that maintain hPSCs in pluripotency (Boyer et al., 2005). We previously showed that PAX6 is necessary and sufficient to generate a human nEc cell fate (Zhang et al., 2010;Chen et al., 2018). Both T (Brachury) and MIXL1 TFs are tightly associated with Me, whereas GATA4 and GATA6 are essential for the initiation of an En fate (Tam and Loebel, 2007;Faial et al., 2015).

Genetic engineering of lineage reporter for high-throughput screening in hPSCs
Within the three germ layers, nEc specification is independent of ectopic signaling pathways activation, referred to as a ''default'' cell fate induction from pluripotent stem cells (Munoz-Sanjuan and Brivanlou, 2002). We have revealed that the PAX6 TF is a human nEc cell fate determinant and thus can serve as a human nEc hallmark gene (Zhang et al., 2010;Chen et al., 2018). Moreover, nEc can be generated directly from hPSCs via overexpression of PAX6, even under strict hPSC culture conditions (Zhang et al., 2010). During early embryonic development, the expression of PAX6 is intricately regulated with cis-regulatory DNA machineries encompassing large genomic regions (Kleinjan et al., 2006;Tyas et al., 2006). We, therefore, engineered hESCs with a PAX6 reporter via CRISPR/Cas9-mediated homologous recombination (HR) to indicate PAX6 expression profiles and screen for intrinsic pathways that either drive hPSCs to three germ layer differentiations or to a nEc fate without changing culture conditions. We first constructed a donor plasmid harboring an in-frame coding cassette of Neo-P2A-tdTomato-T2A-HSVtk flanked by 5 0 AND-3 0 homology arms of PAX6 ( Figure 1A). A single guide RNA (gRNA) targeting exon 4 of PAX6 DNA near the ATG start codon was designed, and the cleavage efficacy was verified in HEK293 cells after co-transfection with Cas9 expression plasmids ( Figure S1A). Sanger sequencing results showed that all 5 predicted loci which showed high similarity with the target sequences were identified intact with no mutations, suggesting no off-target effects of the designed gRNA (Figures S1B and S1C). After correctly HR, the putative integrated hESCs would express Neo-P2A-tdTomato-T2A-HSVtk protein under the tight control of endogenous PAX6 gene expression regulatory machineries. After protease cleavage, tdTomato mirrors an endogenous PAX6 expression and nEc fate, whereas neomycin (Neo) and HSVtk could serve as positive selection and negative depletion tools, respectively, when neomycin or ganciclovir are supplied to the cells. Genomic polymerase chain reaction (PCR) revealed that clones 2 and 11 were two monoallelic HR lines ( Figures 1A and 1B). Southern blot analysis confirmed precise HR with no off-target recombination ( Figure 1C). Both clone 2 and clone 11 retained typical hESC morphology and growth behavior compared with wild-type hESCs. After guiding the lines to a nEc fate with our standard differentiation protocol (Zhang and Zhang, 2010;Chen et al., 2018;Zhang et al., 2001;Li et al., 2009;Liu et al., 2019;Chi et al., 2016), both lines were normally specified to columnar neuroepithelial cells arranged in neural rosette-like structures, and cells started to express tdTomato from day 6 ( Figure S1D). Immunostaining further revealed that tdTomato was co-labeled with endogenous PAX6 in the differentiated nEc cells ( Figure 1D). Therefore, we successfully generated a nEc reporter line in hPSCs (referred to herein as PAX6-tdTomato reporter line), which can be used for lineage specification studies. We modified the two-component inducible Cas9 (iCas9) system (Gonzalez et al., 2014) in the PAX6-tdTomato line. One donor plasmid containing a doxycycline-inducible Cas9 expression cassette (33TRE-FLAG-Cas9) and another carrying a constitutive reverse tetracycline transactivator (M2rtTA) expression cassette (CAG-M2rtTA) were constructed ( Figure 2A). Both donor plasmids were electroporated into the clone 11 PAX6-tdTomato reporter line together with left and right TALEN plasmids targeting AAVS1 (also known as PPP1R12C) loci. After genomic DNA PCR, we identified that clone 4 had both transgenes integrated in the AASV1 loci biallelically ( Figure 2B). Western blotting and immunostaining confirmed the robust induction of Cas9 expression upon doxycycline treatment ( Figures 2C and 2D). In addition, the constructed PAX6-tdTomato/iCas9 line showed typical hPSCs morphology and OCT4 and NANOG expression, as expected ( Figures 2D and S2A). After doxycycline treatment, over 90% cells showed both NANOG and Cas9 staining, suggesting high induction efficiency ( Figures S2B and S2C). (B) Genomic DNA PCR results of individual colonies after gene targeting. Primer sets of PAX6-F1 and PAX6-R1 were used for detecting the wild-type allele, and primer sets of PAX6-F1 and Neo-R were used for detecting the HR allele. Colonies #2 and #11 were monoallelically targeted. (C) Southern blot analysis showed precise but no off-target recombination.

OPEN ACCESS
iScience 24, 101926, January 22, 2021 3 iScience Article To determine the capacity of the iCas9 hPSCs for DNA cleavage, we designed gRNA pairs targeting the human NF1 gene (Liu et al., 2016). The PAX6-tdTomato/iCas9 hESCs were then electroporated with the double gRNAs of NF1. After treatment with doxycycline, we observed 70% (5/7) knockout efficiency based on genomic DNA PCR ( Figure S2D). We also designed a gRNA targeting P53 gene (Liu et al., 2016). The gRNA was packaged into lentivirus, which was then used to infect the PAX6-tdTomato/iCas9 cells. We observed prominent occurrence of indels surrounding the gRNA-guided cleavage site when doxycycline was supplied ( Figure S2E). These results indicate that the iCas9 system constructed through HR at the genetically open AAVS1 loci is tight and functional and thus is suitable for large-scale CRISPR screening.  We hypothesized that genetic loss of function in the PAX6-tdTomato line which caused tdTomato expression would indicate the role of this gene in preventing either trilineage or nEc specification. Therefore, we performed genome-wide loss-of-function screening using the high-throughput CRISPR/Cas9 platform in the PAX6-tdTomato/iCas9 line ( Figure 3A). Lentiviral libraries expressing gRNAs (GeCKO) targeting all coding genes and microRNAs were prepared (Shalem et al., 2014). To avoid spontaneous differentiation, the PAX6-tdTomato/iCas9 hESCs were digested into small clusters (not single cells) and infected with lentiviral libraries at a multiplicity of infection (MOI) of 0.3 to ensure high enrichment in the proportion of cells infected with only one viral particle and thus target a unique gene. As the gRNA lentiviral libraries had a puromycin drug-resistant gene, we removed non-infected cells by 4 days puromycin treatment. The . Genome-scale CRISPR screening in PAX6-tdTomato/iCas9 hESCs (A) Schematic illustration of the genome-scale CRISPR screening strategy. PAX6-tdTomato/iCas9 hESCs were transduced with lentiviral particles of either A pool or B pool gRNA libraries at an MOI of 0.3. Puromycin was applied for 4 days to eliminate non-infected cells. Doxycycline treatment for another 5 days led to target gene loss of functions, and loss of function of germ layer specification preventing genes therefore caused lineage differentiation and aberrant expression of tdTomato under hPSC culture conditions. Wild-type (WT) cells and PAX6-tdTomato cells with no doxycycline treatment did not show tdTomato expression, and they served as negative controls. tdTomato-positive cells induced by doxycycline were purified by FACS. DNA isolation, two-step PCR amplification, and next-generation sequencing were subsequently used for elucidating gRNA enrichment in tdTomato-positive cells.
(B) Epithelial cell clusters expressing tdTomato within the hESCs colony appeared after doxycycline treatment. Scale bar, 100 mm.
(C) Counts of colonies comprising at least one tdTomato-positive cluster of a 6-well plate with or without doxycycline treatment.
(D) FACS of tdTomato-positive cells after doxycycline treatment.
(E) Fluorescent images showed that FACS-sorted cells had uniform tdTomato expression. Scale bar, 100 mm.

OPEN ACCESS
iScience 24, 101926, January 22, 2021 5 iScience Article Figure 4. Identification of lineage-specification prevention genes in hPSCs (A) MAGeCK algorithm was used to estimate the statistical significance of enriched gRNA in tdTomato-positive cells after doxycycline induction. gRNA-associated genes were ranked through robust rank aggregation (RRA) and were termed lineage preventing genes (LPGs). Log2 fold change (LFC) was used to calculate the enrichment index as compared with the library. Each dot of the scatterplot represented one single LPG and genes known for regulating self-renewal of hPSCs were labeled in red. (B) Over 90% of LPGs were efficiently enriched in the second sub-pool library screening.
(C) Percentage of overlapped genes of LPGs and reported essential genes (EGs) related to hPSC proliferation and survival or hPSC self-renewal-required genes. Gene Ontology (GO) analyses showed these overlapped genes were mostly related to mRNA processing and transcriptional regulation of pluripotency, with representative genes of each GO term been boxed and marked in red. iScience Article PAX6-tdTomato/iCas9 hESCs infected with the gRNA lentiviral libraries were subsequently treated with doxycycline for another 5 days. Without doxycycline treatment, we did not see tdTomato-positive cells containing colonies, suggesting high quality hESC maintenance under strict culture conditions. Furthermore, doxycycline treatment repetitively yielded tdTomato-positive clusters within the hESC colonies at a frequency of 4%-5% ( Figures 3B and 3C). Notably, the tdTomato-expressing cells showed columnar neuroepithelial morphology, forming clusters of similar size and containing distinct boarders with surrounding tdTomato-negative hPSCs ( Figure 3B). These tdTomato-positive cells completely lost pluripotent marker OCT4 but were positive for SOX2, a known TF also for facilitating nEc fate decision (Wang et al., 2012) (Figure S3A). These results indicate that the tdTomato-expressed clusters were guided to a nEc fate and their uniform clonal appearance implied a single-cell origin likely caused by a loss-of-function mutation of a gene. Fluorescence-activated cell sorting (FACS) experiment showed that approximately 0.5% of entire hESC population was tdTomato positive after doxycycline treatment ( Figure 3D). We then enriched the tdTomato-positive cells of interest using FACS and validated their tdTomato expression under the microscope immediately and their nEc characteristics by analyzing the expression of nEc marker genes including PAX6, SOX1, DLK1, MEIS2, and ZIC1 at the mRNA level ( Figures 3E and S3B). After validation, the total genomic DNA of these tdTomato-positive cells was extracted, the fragments containing the integrated gRNA were amplified via PCR, and target gRNAs were profiled through sequencing.

Recovered gRNAs are enriched in biological processes essential for hPSC maintenance and lineage specification
To evaluate the CRISPR/Cas9 screening performance, we first verified the integrality of the gRNA pools (A pool and B pool as designed) transduced in the PAX6-tdTomato/iCas9 hESCs without doxycycline treatment, with results showing more than 90% coverage ( Figure S4A). The Gini index was then calculated to reflect disparities in the distribution of gRNAs. Both the A and B pool in the doxycycline-untreated cells showed a low Gini index value, thus representing highly uniform transduction with no bias ( Figure S4B). In contrast, gRNAs retrieved from doxycycline-induced tdTomato-positive cells showed a high Gini index value, indicating target gene selectivity and enrichment ( Figure S4B). Furthermore, high repeatability between independent screens reinforced the reliability of the forward screen system based on lineage reporter and inducible gene knockout ( Figure S4C).
The MAGeCK algorithm (Li et al., 2014) was used to score the statistical significance (using a negative binomial test) of enrichment for individual gRNAs in tdTomato-positive cells by comparison with the overall gRNA pool. Significantly, enriched gRNAs were defined with an over 2-fold change and false discovery rate (FDR) below 0.05 ( Figure S5A). Genes considered relevant were identified by their targeting gRNAs shown in the screen ranked consistently higher (by significance) using robust rank aggregation (RRA). Based on these criteria, we identified 4% (829/19,062) of overall genes in the coding genome as essential and their loss of function in hPSCs triggered lineage specification ( Figure 4A). We therefore termed these genes as lineage-specification preventing genes (LPGs). To further validate these LPGs, we designed a subpool library, which included 182 non-targeting control gRNAs, and 7,324 gRNAs for the 829 LPGs and for those genes which were identified for at least one screen but showed no statistical significance. The subpool gRNAs were then verified by next-generation sequencing and showed a Gini index value at 0.029, suggesting uniform distribution ( Figure S5B). Two sub-pool validation screenings were then performed, and more than 90% of previously identified LPGs were preferentially enriched with over 2-fold changes and FDRs below 0.05 ( Figure 4B). The top 100 LPGs were listed in Table S1.
Both OCT4 and NANOG were among the LPGs, coinciding with their pivotal role in pluripotency (Boyer et al., 2005). TGFBR1, receptor of the Activin/Nodal signaling pathway, which is essential for hPSC maintenance, was also among the identified LPGs (Vallier et al., 2005;Xu et al., 2008). In addition, HDAC3 and EP300, which are known to facilitate self-renewal while inhibiting lineage specification of hPSCs, were also detected within the LPGs (Qiao et al., 2015). Importantly, the identified LPGs were distributed across all chromosomes without enrichment in specific chromosomal regions ( Figure S5C). These results suggest Signaling molecules related to Wnt, TGF-b, FGF/MAPK, and Notch pathways are early embryonic developmental cues, which showed a tight link with pluripotent genes and lineage specifiers. This PPI analysis pointed to an antagonistic effect within early lineage specifiers in hPSCs, given loss of function of mesendoderm specifiers led to spontaneous nEc differentiation.

OPEN ACCESS
iScience 24, 101926, January 22, 2021 7 iScience Article that the genome-wide CRISPR screening tool revealed a plethora of candidate LPGs and the mutation of which caused hPSC pluripotency exit and germ layer entry.
EGs critical for pluripotent stem cell survival and proliferation have been defined in both diploid and haploid hPSCs through CRISPR loss-of-function screening (Yilmaz et al., 2018;Mair et al., 2019). We made an intersection between LPGs and EGs and found that less than 10% of LPGs overlapped with the EGs, and these overlapped genes were enriched in the biological function of RNA processing ( Figure 4C). With the OCT4-GFP reporter, genes required for hPSC self-renewal were also identified through RNAi library screening (Chia et al., 2010). LPGs shared far more genes with self-renewal-required genes than with the common genes of LPGs and EGs ( Figures 4C and S5D). In addition, the shared LPGs and self-renewal-required genes were mostly clustered in the transcriptional regulation of pluripotency ( Figures 4C and S5D). These results suggested that our CRISPR screening successfully retrieved a group of genes essential for hPSC self-renewal. However, given the advanced lineage reporter design and forward screening system, we also discovered those previously undefined genes essential for other biological processes, which shape trilineage specification.
We further annotated the LPGs to determine their biological function using the Panther Classification System (Mi et al., 2005). Biological process analysis showed strong enrichment in metabolic processes, biological regulation, and cellular component organization pathways ( Figure S5E). Protein classification analysis showed enrichment in nucleic acid binding, transcription factors, enzyme modulators, and signaling molecules ( Figure S5F). Additionally, we assigned various functional categories to hierarchical clusters based on Kappa-statistical similarities among their gene memberships using Metascape (Tripathi et al., 2015). This revealed that the LPGs were mostly categorized into eight functional modules, among which embryonic development, epigenetic regulation, metabolic regulation, and mRNA processing were the most prominent ( Figures 4D and S5G).
To further delineate the regulatory networks controlling LPG expression, we profiled the DNase hypersensitive sites (DHSs) for these genes. DHSs are regions of open chromatin that make DNA accessible, and these accessible chromatin zones are functionally related to transcriptional activity via recruitment of dominant TFs. Here, compared with random or background genes, LPGs harbored significantly higher DHS coverage through chi-square tests ( Figures 4E and 4F), suggesting their open chromatin states and active or poised transcription. The DNA sequences around the DHSs of all LPGs were subsequently subjected to TF binding motif analysis. We found that the top motifs were well-characterized TF binding sites associated with pluripotency, such as SP1 of the Z-finger domain, OCT4/NANOG/SOX2 core pluripotency factors, PBX3 of homeobox domain, NF-kB repressing factor (NRF), and Yin Yang-1 (YY1) ( Figure 4G). These results suggest that LPGs are largely transcriptionally active genes in hPSCs and are tightly regulated by core pluripotency factors, reinforcing the fidelity and biological relevance of the identified LPGs. Although we found specific motif enrichment in modules such as metabolism, autophagy, and immune response, the SP1, OCT4/NANOG/SOX2, and YY1 binding motifs appeared in multiple modules ( Figure S5H), reflecting their broad effect in regulating pluripotency and lineage specification.

Mesendodermal lineage specifiers antagonize nEc cell fate
Those LPGs within the embryonic development module could be clearly separated into a large-sized and a small-sized group. The large-sized group contained core pluripotency TFs, such as OCT4 and NANOG, and key molecules related to TGF-b, FGF/MAPK, and Notch signaling, which largely overlapped with pluripotency genes retrieved from previous RNAi screening using the OCT4-GFP reporter (Chia et al., 2010). Strikingly, the small-sized group was centered by typical lineage specifiers, such as T and MIXL1 for Me and GATA4 and EOMES for En. Protein-protein interaction (PPI) network analysis showed the interplay of these pluripotency genes and Me or En specifiers ( Figure 4H). Deng and colleagues revealed that a combination of nEc, Me, and En lineage specifiers can efficiently reprogram fibroblasts into iPSCs (Shu et al., 2013). It has been reported that germ layer genes and pluripotent genes are co-expressed in a small population of hPSCs, and these cells still preserved their clonogenicity and potency to differentiate into all three germ layers (Allison et al., 2018;Han et al., 2010;Hough et al., 2014). Thus, expression of trilineage specifiers in hPSCs could act against each other to maintain a state of equilibrium in hPSCs, and upsetting such a balance will cause the differentiation of hPSCs toward antagonistic lineage. iScience Article epigenetic regulation in LPGs, mostly categorized into families of histone acetylation or de-acetylation, methylation or de-methylation, chromatin remodelers, and others related to ubiquitination or phosphorylation ( Figure 5A). Many of these genes were downregulated during lineage specification ( Figure S6A). To understand the roles of these epigenetic regulations in trilineage specification, we profiled the histone modifications (HMs) in gene categories of hPSCs-, mesendoderm-, and neural lineage-specific genes. (D) miRNAs within the scope of LPGs were identified by MAGeCK algorithm. The significance of each miRNA was ranked through RRA. Each dot of the scatterplot represented one single miRNA and top 5 miRNA, as well as hsa-mir-302e/b known for regulating self-renewal of hPSCs were labeled in red.

Mutations of epigenetic machineries lead to hPSC pluripotency exit or lineage differentiation
(E) Among all miRNAs in miRBae, hPSCs-, Ec-, or Me&En-related miRNAs were defined by using PCA classifier according to the enrichment score of all 3 categories based on GO annotation of its target genes. Red dots represented miRNAs shown in defined LPGs, and they were highly enriched in the Ec-related group.
(F) Motif analyses with HOMER for the promoter regions of Ec-or Me&En-related miRNAs.

OPEN ACCESS
iScience 24, 101926, January 22, 2021 9 iScience Article Results showed that during mesendoderm and nEc lineage specification, the levels of histone acetylation/ H3K4me3 in promoter regions and H3K4me1/H3K27ac in enhancer regions of lineage-specific genes were significantly increased, whereas those of H3K27me3 were decreased ( Figures 5B and S6B). Heatmap analysis confirmed an obvious decrease in histone deacetylases, such as BRMS1, MTA3, HDAC3, and SKI, and methyl transferases for H3K27, including EHMT1 and MTF2 along with trilineage differentiation of hPSCs ( Figure 5C). These results suggest that active epigenetic machineries are repressed, while repressive epigenetic machineries are highly active in lineage genes in hPSCs. Thus, loss of function of genes centered in these epigenetic machineries caused aberrant activation of lineage genes in hPSCs and therefore led to spontaneous lineage differentiation ( Figure S6C).
hPSCs-expressed miRNAs are clustered and maintain pluripotency by targeting trilineagespecific genes We identified a cluster of miRNAs (85/1864) within the screened LPGs ( Figure 5D). Among these miRNAs, the hsa-mir-302 family is reported to play a role in hPSC maintenance by targeting COUP-TFII, a TF expressed in the nEc (Rosa and Brivanlou, 2011). To profile the function of miRNAs in LPGs, we annotated the top five miRNAs through miRBase (Griffiths-Jones et al., 2006) and performed Gene Ontology (GO) analysis for their target genes. We found that most of the target genes were highly enriched in the nEc lineage and were associated with nEc development ( Figures S6D and S6E). Indeed, when we sub-grouped the miRNA superfamily into hPSC-, mesendoderm-, and ectoderm-related groups according to their target genes, we found that most miRNAs in our defined LPGs were in the Ec-related group ( Figure 5E). These data suggested that hPSC-expressed miRNAs act as differentiation barriers, which specifically inhibit differentiation of an individual lineage. In line with this hypothesis, the promoter regions of Ec-related miRNAs primarily exhibited binding sites for Zinc finger family TFs and NANOG, whereas the mesendoderm-related miRNAs were enriched with distinct binding motifs of TFs, such as SOX2/4 and forkhead family TFs ( Figure 5F).

Metabolic pathways profoundly regulate hPSC maintenance and lineage specification
An accumulation of evidence has demonstrated that metabolic pathways play important roles in hPSC maintenance (Wu et al., 2016). Here, genes belonging to the metabolic pathways were the largest superfamily in the identified LPGs. GO analysis showed that these metabolism-related genes were mostly clustered in functional modules of the TCA cycle, energy production, and lipid or amino acid metabolism ( Figure 6A).
While hPSCs generally use glycolysis for energy production, a switch from glycolysis to OXPHOS is observed during their differentiation (Zhang et al., 2016). Within the LPGs, we were surprised to find a cluster of mitochondrial electron transport complex genes, considered as the endpoint of the TCA cycle. More importantly, these mitochondrial TCA-related genes were highly expressed in undifferentiated hPSCs but exhibited decreased expression upon differentiation ( Figure 6B), implying their important role in hPSCs. They may act like c-Myc, an important modulator for mitochondrial biogenesis and is key for hPSC pluripotency maintenance and hiPSC reprogramming (Takahashi et al., 2007;Cliff et al., 2017). Although systematic research is still required to fully elucidate the impact of bioenergetic pathways on the fate of hPSCs, our data suggest that the mitochondrial TCA cycle and OXPHOS are absolutely required for the self-renewal of hPSCs, as loss of function of genes in these pathways caused spontaneous lineage differentiation.
We also found that most LPGs related to amino acid metabolism were downregulated during the process of hPSC differentiation ( Figure 6C), implying an important role of amino acid metabolism in maintaining hPSCs (Shiraki et al., 2014). Surprisingly, genes related to cholesterol synthesis were highly enriched in LPGs (8 genes). These cholesterol synthesis genes were highly expressed in undifferentiated hPSCs and Me or En lineages but were downregulated when a nEc cell fate was adopted ( Figure 6D). This implies that, different from other metabolic pathways, the cholesterol synthesis pathway is a preset barrier and maintains hPSCs via retardation of nEc specification.

Discussion
With a nEc-specific reporter and inducible CRISPR/Cas9 expression system, we performed genome-scale loss-of-function screening in hPSCs to identify dominant genes that determine pluripotency and lineage ll OPEN ACCESS 10 iScience 24, 101926, January 22, 2021 iScience Article Figure 6. Cholesterol synthesis pathway restricts nEc specification from hPSCs (A) PPI network analysis generated by STRING revealed that LPGs within the metabolism module comprised key molecules function in multiple metabolic pathways. (B) Genes related to TCA cycle and oxidative phosphorylation were enriched in LPGs, and heatmap showed decreased expression of these LPGs along trilineage differentiation of hPSCs. The expression level was normalized to RPKM and z-scored.
(C) Heatmap showed decreased expression of LPGs related to amino acid metabolism along trilineage differentiation of hPSCs. The expression level was normalized to RPKM and z-scored.
(D) Genes of cholesterol synthesis pathway that were enriched in LPGs were labeled in red. Heatmap showed that the expression of these cholesterol synthesis pathway-related LPGs was specifically downregulated when the hPSCs were specified toward a nEc fate. The expression level was normalized to RPKM and z-scored.
(E) A schematic model for hPSC self-renewal and trilineage differentiation. Extracellular signaling, epigenetic regulation, transcriptional regulation, and metabolic pathways play a synergistic role in regulating hPSC self-renewal and trilineage differentiation. Histone modifications share common mechanisms in leveraging the transcriptional activity of genes belonging to all three lineages. The TCA cycle, the respiratory chain, and the amino acid metabolism are all important for maintaining an undifferentiated state of hPSCs, though glycolytic flux is highly active. A mutual antagonistic effect occurs in Me&En and nEc lineage specifiers, and a balance of these lineage specifiers is required for the self-renewal of hPSCs. hPSC-enriched miRNAs showed clear lineage gene targeting preference. The cholesterol synthesis pathway plays a pivotal role in maintaining hPSCs via restraining the nEc differentiation.

OPEN ACCESS
iScience 24, 101926, January 22, 2021 11 iScience Article specification. We discovered various LPGs for hPSCs and loss of function of which led to either pluripotency exit or specific lineage commitment. Our study demonstrated that lineage-based reporters combined with CRISPR/Cas9 whole-genome screening offer a powerful strategy for revealing the biological processes and gene networks underlying specific cell fate determination, which is a key question in developmental biology and the basis for regenerative medicine. The identified LPGs could also serve as an invaluable resource for future studies on cellular and molecular events regarding hPSC self-renewal and trilineage development.
Our study re-emphasized the pivotal role of Notch/TGF-b/FGF signaling and OCT4/NANOG core transcriptional networks in regulating pluripotency of hPSCs (Vallier et al., 2005;Xu et al., 2008;Boyer et al., 2005). Our screening also revealed that common cellular biological processes, such as cellular organization, cell cycle regulation, and RNA processing, also play universal and important roles in hPSCs. Some of these functional modules overlap with EGs crucial for hPSC survival and proliferation. As demonstrated by our screening, there was a specific module for immunological responses within the LPGs. Eggenberger et al. showed that ectopic activation of the canonical type I interferon antiviral pathway drives hPSCs away from pluripotency and induces aberrant trilineage differentiation, supporting the conclusion that immune responses are also tightly associated with pluripotency and lineage development (Eggenberger et al., 2019).
Epigenetic modifiers establish a fundamental genomic environment for precise gene expression control and cellular context homeostasis. We identified that HMs and chromatin remodelers had profound effects on hPSC self-renewal and lineage differentiation, based on the following observations: (1) HM-related machineries and chromatin remodelers in LPGs were largely downregulated during lineage specification; (2) active HMs in the promoter and enhancer regions of lineage specific genes were enriched, whereas repressive HMs were de-enriched during corresponding lineage commitment; thus we concluded that HMs and chromatin remodelers maintained hPSCs in an undifferentiated state by establishing transcriptional barriers to restrict lineage gene expression. Notably, those miRNAs highly expressed in hPSCs retained pluripotency by targeting lineage genes. However, these miRNAs could be clearly sub-grouped into specific lineage-related clusters. Moreover, the promoter regions of the ectoderm-related miRNAs were enriched in binding sites of NANOG, the core pluripotency factor that antagonizes Ec but favors definitive En development. With our analysis, we also found that the mesendoderm-related miRNAs were enriched with binding motifs for SOX2, a known TF for facilitating nEc but inhibiting Me or En specification (Wang et al., 2012).
A small but intriguing cluster within the identified LPGs was the key mesendoderm specifiers, i.e., GATA4, EOMES, MIXL1, and T. Deng and his colleagues revealed that GATA3 and DLX3 can substitute OCT4 and SOX2 to reprogram mouse fibroblasts into iPSCs. Furthermore, they established a ''seesaw'' model, which demonstrated the balance between pluripotency factors and/or counteracting lineage specifiers (Shu et al., 2013). The same model also applies in humans, with lineage-specific genes able to efficiently reprogram human fibroblasts into hiPSCs (Mai et al., 2018;Shu et al., 2015). With the advent of single-cell RNAseq technology, pluripotent stem cells are considered relatively heterogeneous, and lineage-specific genes are also expressed in some subpopulations of pluripotent stem cells responsible for clonal expansion (Allison et al., 2018;Hough et al., 2014). Thus, we hypothesized that, at least, in some populations of pluripotent stem cells, lineage specifiers counteract with each other and therefore maintain their pluripotency. In the future, the same strategy by applying mesendoderm-specific reporters could be used to identify Ec including nEc specifiers, LPGs antagonizing mesendoderm differentiation.
Here, genes within the metabolic pathway accounted for the largest family among identified LPGs. Although glycolytic flux has been shown to be crucial for the maintenance of hPSCs (Folmes et al., 2011;Panopoulos et al., 2012), we identified LPGs that were clustered in the functional module of mitochondrial TCA cycle and OXPHOS. Multiple studies have shown that proteins function in the TCA cycle and OXPHOS are higher in hPSCs than in their differentiated derivatives (Varum et al., 2011). Thus, the mitochondrial TCA cycle and OXPHOS and are also indispensable for hPSC pluripotency. Surprisingly, loss of function of genes within the cholesterol synthesis pathway led to nEc differentiation of hPSCs. It has been reported that the expression of pluripotent gene OCT4 was decreased when the cholesterol biosynthesis pathway was inhibited (Matsuzaki et al., 2018). Moreover, the cholesterol synthesis pathway is highly active in hPSCs and Me and En lineages but low in nEc. Thus, both glycolysis and OXPHOS bioenergetic pathways are crucial for the maintenance of hPSCs, and cholesterol synthesis process is required for human pluripotency via retardation of nEc specification. In summary, this study provides a full landscape of biological processes underlying hPSC self-renewal and trilineage differentiation ( Figure 6E). Pathways involved in bioenergetic production, HMs and chromatin remodelers, immune responses, RNA processing, as well as core pluripotency-related signaling and transcriptional networks, function synergistically to maintain a pluripotent state in hPSCs via inhibition of trilineage differentiation. However, hPSC-enriched miRNAs, lineage specifiers, and the cholesterol synthesis pathway exhibit germ layer selectivity, and these biological processes promote or antagonize specific germ layers along the maintenance of hPSCs or early lineage differentiation.

Limitations of the study
In this paper, we performed CRISPR/Cas9 screen in hPSCs to identify those genes that prevent lineage specification. However, when taking the trilineage cell fate decision in count, we only used the PAX6 reporter and lacked parallel screen with the mesendodermal reporters. The uncovered library of lineage specification preventing genes is a valuable source for studying mechanisms underlying hPSC self-renewal and trilineage specification. However, specific validation is apparently required when comes to individual genes.

Materials availability
Requests for materials and reagents should be directed to the Lead contact.

Data availability
All relevant data are available from the correspondence authors upon reasonable request. The accession number for the CRISPR/Cas9 screening data reported in this paper is GEO: GSE132309.

Methods
All methods can be found in the accompanying Transparent Methods supplemental file.

Declaration of interests
The authors declare no competing interests.        shown in LPGs demonstrated that they were preferentially nEc development related.

Neural Differentiation
Detailed procedure for embryoid body (EB)-based differentiation of hPSCs toward a PAX6 positive nEc and cortical neural progenitors was described previously (Zhang and Zhang, 2010;Chi et al., 2016;Zhang et al., 2001). Briefly, hESCs were digested with dispase, lifted from the MEF layers, and triturated into 100-200 μm pieces. EBs formed after 4 days of suspension culture in hESCM.
EBs were then transferred into the neural induction medium (NIM) (DMEM/F12 with 1×NEAA, 1×N2 supplement and 2 mg/ml heparin) for another two days and then plated on laminin-coated culture surfaces in NIM. nEc formed at day 6-10 post differentiation.

Lentivirus Packaging
HEK 293FT cells (Invitrogen) of 70% to 80% confluence in a 10 cm dish were transiently transfected with 10 μg gRNA lentiviral vectors, 5 μg VSVG and 7.5 μg PAX2 via the calcium phosphate precipitation method. Fresh DMEM with 10% FBS was then supplied 14 hr after transfection. 48 hr later, lentivirus particles-containing supernatant was harvested and filtered through a 0.45 μm cellulose acetate filter, followed by ultracentrifugation (Beckman) with a speed of 55000 g at 16 °C for 3 hours. Lentivirus particles-containing pellet was reconstituted in hESCM.
The gRNA GCCCCATATTCGAGCCCCG sequence, targeting a region close to the ATG start codon of human PAX6 gene was designed through the Zhanglab website (https://zlab.bio/guidedesign-resources) and cloned into the gRNA expressing vector as previously described (Chen et al., 2018;Chi et al., 2016). Cleavage efficiency of the designed gRNA was tested in HEK 293FT cells through co-transfection with Cas9 expressing plasmids (Addgene #44719). After 60~65 hr of transient transfection, genomic DNA was extracted and the DNA fragments surrounding the cleavage site of PAX6 gene was PCR amplified for Sanger sequencing. Overlapped peaks indicate efficient double strand breaks made by the cleavage followed by indels generated after DNA repair.

Genome-wide CRISPR Screening
The lentiviral gRNA library for genome-wide CRISPR/Cas9 screening was purchased from were then gel purified, quantified, and sequenced using the pair-ended 150 bp sequencing protocol (Illumina Hiseq-x-ten system).

Sub-pool library validation screen
In order to construct a sub-pool gRNA library for validation, we modified the lentiCRISPR v2 (Addgene, #52961) by replacing Cas9 with blastcidin. 7,506 sgRNAs, including 7,324 selected hit gRNAs from the primary screen targeting 3,941 genes, and 182 non-targeting control sgRNAs, were individually cloned into the modified lentiCRISPR v2 vector and verified by next generation sequencing.
Lentivirus was produced with the sub-pool library. PAX6-tdTomato/iCas9 hESCs were treated with dispase and triturated into small clumps. A total of 1.

Southern Blot Analysis
The DNA probe targeting tdTomato sequence was PCR amplified from the PAX6-tdTomato donor vector with the PCR DIG Probe Synthesis Kit (Roche). Primer sets were as follows: 5'-TGGTGAGCAAGGGCGAGGAG-3' and 5'-TGCCGGTGCTGCCGGTGCCAT-3'. The efficacy of the DIG-labeled probe was evaluated with dot blot. 10 μg of genomic DNA extracted from PAX6-tdTomato hESCs was digested over-night with BsaI and subjected to agarose gel electrophoresis in TBE buffer. Gels were then denatured, neutralized, and transferred over night by capillarity on Hybond-N membranes (GE Healthcare) using 20×SSC transfer buffer. Hybridization was carried out overnight at 50°C and then was blocked with Roche blocking buffer. DIG signals were detected with an AP-conjugated DIG-Antibody (Roche) using CDP-Star (Roche) as a substrate for chemiluminescence.

Data Analysis for CRISPR/iCas9 Screening
We sequenced in pair-ended mode with 150-bp read length, first in read paired ( The count table for each gRNA was then normalized relative to the total number of reads in each of conditions and two libraries were combined. As a result, we got the normalized counts of each gRNA in control sample (unscreened cell library) and KO-screened sample respectively. The MAGeCK algorithm (Li et al., 2014) was used to estimate the statistical significance (using a negative binomial test) of enrichment for each gRNA in the KO-screened group compared to control group. LPGs were then identified by looking for genes whose gRNAs were ranked consistently higher (by significance) using robust rank aggregation (RRA). The negative controls were incorporated in the MAGeCk analysis to generate null distributions and calculate the p-value and FDR for each gene.
The gRNA sequences were aligned to the gRNA sub-pool Burrows-Wheeler index with Bowtie2 (version 2.3.4.3) and the reads with unique alignment were used for the gRNAs reads counting.
Each gRNA were calculated as LFCs of gRNA counts between control sample (unscreened cell library) and KO-screened sample. The gRNAs counts were normalized by comparing to median of non-targeting control gRNAs. We averaged the LFCs of two replicates and genes with an over 2fold change and FDR below 0.05 were positive.

High-throughput Sequencing (NGS) Data Analysis
Data of DNase-seq (DNase I hypersensitive sites sequencing) used to identify the location of regulatory regions was from Gene Expression Omnibus (GEO) database (GSM878612, GSM878613). Regulatory motifs of the DHs of LPGs were predicted by using HOMER software (Heinz et al., 2010). To depict the epigenetic regulation of the committed nEc fate from hESCs, we performed the profiles of various histone modifications by using computeMatrix and plotProfile commands of deepTools (Ramirez et al., 2014), and the ChIP-seq data was from Roadmap Epigenomics project. To reveal the changes of mRNA expression levels of LPGs during different lineage specification from hESCs, we download RNA-seq data form GEO database (Gifford et al., 2013). Expression levels for genes were normalized to Reads Per Kilobase per Million mapped reads (RPKM) and row z-scored when performing heatmap profiles.

Functional Annotation of Defined LPGs
We performed functional classification of LPGs by using the Panther Classification System (Mi et al., 2005). And then we assigned various functional categories to clusters of genes using Metascape (Tripathi et al., 2015) with default analysis parameters. In details, all the statistically enriched terms identified from Gene Ontology (GO), KEGG pathway and Reactome Gene Sets were hierarchically clustered into a tree based on Kappa-statistical similarities among their gene memberships. Clusters which represent a group of similar functional categories were defined by a kappa score of 0.3 as threshold. A subset of representative terms from each cluster was automatically selected by Metascape and converted into a network, where terms with similarity above 0.3 are connected by edges. More specifically, terms with the best p-values from each of the clusters were depicted as network nodes, with the constraint of having a maximum of 15 terms per cluster and 250 terms in total.

Gene Network Construction
The online database resource search tool (STRING) for retrieving interactive genes provides protein-protein interaction information, including the prediction and comparison of inter-genomic interactions (Franceschini et al., 2013). Only protein associations with combined confidence score above 0.5 were retained and used to assign weights to each network link. Different clusters of geneencoded proteins were represented by nodes of different colors in the network, and line shape indicated the predicted mode of action between proteins referring to https://string-db.org/ (version