Search

Measuring chromatin accessibility is a powerful method to identify cell types and states. Performed at single-cell resolution, this assay has generated catalogs of genome-wide DNA regulatory sites, whole-organism cell atlases, and dynamic chromatin reorganization through development. However, the limited throughput of current single-cell approaches poses a challenge for implementing proper study designs, population-scale profiling, and/or very deep profiling of complex samples. To this end, we developed a 10X-compatible combinatorial indexing ATAC sequencing (“txci-ATAC-seq”), which is a combinatorial indexing framework that initially indexes (“pre-indexes”) chromatin within nuclei with barcoded transposases followed by encapsulation and further barcoding using a commercialized droplet-based microfluidics platform (10X Genomics). Leveraging this molecular hashing strategy, we demonstrate that txci-ATAC-seq enables the indexing of up to 200,000 nuclei across multiple samples in a single emulsion reaction, representing a ∼22-fold increase in throughput compared to the standard workflow at the same collision rate. To improve the efficiency of this new technique, we further developed a faster version of the protocol (“Fast-txci-ATAC-seq”) that separates sample pre-processing from library generation and has the potential to profile up to 96 samples simultaneously. We initially benchmarked our assay by generating chromatin accessibility profiles for 230,018 cells from five native tissues across three experiments, including human cortex (28,513 cells), mouse brain (48,997 cells), human lung (15,799 cells), mouse lung (73,280 cells), and mouse liver (63,429 cells). We also applied our method to a club cell secretory protein knockout (CC16-/-) mouse model to examine the biological and technical limitations of the mouse line. By characterizing DNA regulatory landscapes in 76,498 wild-type and 77,638 CC16-/-murine lung nuclei, our investigations uncovered previously unappreciated residual genetic deviations from the reference strain that resulted from the method of gene targeting, which employed embryonic stem cells from the 129 strain. We found that these genetic remnants from the 129 strain led to profound cell-type-specific changes in chromatin accessibility in regulatory elements near a host of genes. Collectively, we defined single-cell chromatin signatures in 384,154 nuclei from 13 primary samples across different species, organs, biological replicates, and genetic backgrounds, establishing txci-ATAC-seq as a robust, high-quality, and highly multiplexable single-cell assay for large-scale chromatin studies.

Gene regulatory networks, which control gene expression patterns in development and in response to stimuli, use regulatory logic modules to coordinate inputs and outputs. One example of a regulatory logic module is the gene regulatory cascade (GRC), where a series of transcription factor genes turn on in order. Synthetic biologists have derived artificial systems that encode regulatory rules, including GRCs. Furthermore, the development of single-cell approaches has enabled the discovery of gene regulatory modules in a variety of experimental settings. However, the tools available for validating these observations remain limited. Based on a synthetic GRC using DNA cutting-defective Cas9 (dCas9), we designed and implemented an alternative synthetic GRC utilizing DNA cutting-defective Cas12a (dCas12a). Comparing the ability of these two systems to express a fluorescent reporter, the dCas9 system was initially more active, while the dCas12a system was more streamlined. Investigating the influence of individual components of the systems identified nuclear localization as a major driver of differences in activity. Improving nuclear localization for the dCas12a system resulted in 1.5-fold more reporter-positive cells and a 15-fold increase in reporter intensity relative to the dCas9 system. We call this optimized system the “Synthetic Gene Regulatory Network” (SGRN, pronounced “sojourn”).

ObjectiveEndothelial cells (ECs) regulate atherogenesis with Endothelial-to-Mesenchymal (EndMT) transition correlating with disease. Single cell (sc) effects of EndMT perturbationsin vitro,with quantitative comparison to signatures of human lesionsin vivois lacking.Approach and ResultsMultiomic profiling that concurrently measures sc RNA and ATAC was performed on 8 distinct primary cultures of human aortic ECs exposed to 7-day of EndMT-promoting perturbations: IL1B, TGFB2, and ERG knockdown (siERG). Meta-analysis of sc transcriptomes across 17 human arterial specimens was performed and two quantitative measures assessed the similarities ofex vivoversusin vitromolecular profiles. Primary HAEC cultures were reproducibly populated by 4 major clusters, termed EC1 -4: EC1-angiogenic; EC2-proliferative; EC3-activated/mesenchymal-like; and EC4-mesenchymal. Independent exposure to siERG, IL1B and TGFB2 elicited mostly distinct transcriptional and chromatin accessible responses. EC1 and EC2, the most canonically ‘healthy’ EC populations were affected predominantly by siERG; the activated cluster EC3 was most responsive to IL1B; and the mesenchymal population EC4 was most affected by TGFB2. Quantitative comparisons betweenin vitroandex vivotranscriptomes confirmed EC1 and EC2 as most canonically EC-like, and EC4 as most mesenchymal with minimal effects elicited by siERG, IL1B and TGFB2. Lastly, accessible chromatin regions unique to EC2 and EC4 were most enriched for CAD-associated SNPs from GWAS suggesting these cell phenotypes harbor CAD-modulating mechanisms.ConclusionPrimary EC cultures contain markedly heterogenous cell subtypes defined by their molecular profiles. Surprisingly, pro-EndMT exposures for 7 days were inadequate to shift cells from one subpopulation to another suggesting relatively stable molecular phenotypes in culture. Interpretations could be that EndMT acts on a modest number of transcripts or that thein vitrosystems used herein fail to recapitulate the complex EndMT-promoting microenvironment of human atherosclerotic lesions. Recognizing and leveraging heterogeneityin vitroshould improve fidelity of these systems for modelingin vivobiology.

Recently developed single-cell technologies allow researchers to characterize cell states at ever greater resolution and scale. Caenorhabditis elegans is a particularly tractable system for studying development, and recent single-cell RNA-seq studies characterized the gene expression patterns for nearly every cell type in the embryo and at the second larval stage (L2). Gene expression patterns give insight about gene function and into the biochemical state of different cell types; recent advances in other single-cell genomics technologies can now also characterize the regulatory context of the genome that gives rise to these gene expression levels at a single-cell resolution. To explore the regulatory DNA of individual cell types in C. elegans, we collected single-cell chromatin accessibility data using the sci-ATAC-seq assay in L2 larvae to match the available single-cell RNA-seq data set. By using a novel implementation of the latent Dirichlet allocation algorithm, we identify 37 clusters of cells that correspond to different cell types in the worm, providing new maps of putative cell type–specific gene regulatory sites, with promise for better understanding of cellular differentiation and gene regulation.

Glycosylphosphatidylinositol-anchored high density lipoprotein–binding protein 1 (GPIHBP1), the protein that shuttles LPL to the capillary lumen, is essential for plasma triglyceride metabolism. When GPIHBP1 is absent, LPL remains stranded within the interstitial spaces and plasma triglyceride hydrolysis is impaired, resulting in severe hypertriglyceridemia. While the functions of GPIHBP1 in intravascular lipolysis are reasonably well understood, no one has yet identified DNA sequences regulating GPIHBP1 expression. In the current studies, we identified an enhancer element located ∼3.6 kb upstream from exon 1 of mouse Gpihbp1. To examine the importance of the enhancer, we used CRISPR/Cas9 genome editing to create mice lacking the enhancer (Gpihbp1(Enh/Enh)). Removing the enhancer reduced Gpihbp1 expression by >90% in the liver and by ∼50% in heart and brown adipose tissue. The reduced expression of GPIHBP1 was insufficient to prevent LPL from reaching the capillary lumen, and it did not lead to hypertriglyceridemia—even when mice were fed a high-fat diet. Compound heterozygotes (Gpihbp1(Enh/−) mice) displayed further reductions in Gpihbp1 expression and exhibited partial mislocalization of LPL (increased amounts of LPL within the interstitial spaces of the heart), but the plasma triglyceride levels were not perturbed. The enhancer element that we identified represents the first insight into DNA sequences controlling Gpihbp1 expression.

Background The “Assay for Transposase Accessible Chromatin sequencing” (ATAC-seq) is an efficient and easy to implement protocol to measure chromatin accessibility that has been widely used in multiple applications studying gene regulation. While several modifications or variants of the protocol have been published since it was first described, there has not yet been an extensive evaluation of the effects of specific protocol choices head-to-head in a consistent experimental setting. In this study, we tested multiple protocol options for major ATAC-seq components (including three reaction buffers, two reaction temperatures, two enzyme sources, and the use of either native or fixed nuclei) in a well-characterized cell line. With all possible combinations of components, we created 24 experimental conditions with four replicates for each (a total of 96 samples). In addition, we tested the 12 native conditions in a primary sample type (mouse lung tissue) with two different input amounts. Through these extensive comparisons, we were able to observe the effect of different ATAC-seq conditions on data quality and to examine the utility and potential redundancy of various quality metrics. Results In general, native samples yielded more peaks (particularly at loci not overlapping transcription start sites) than fixed samples, and the temperature at which the enzymatic reaction was carried out had a major impact on data quality metrics for both fixed and native nuclei. However, the effect of various conditions tested was not always consistent between the native and fixed samples. For example, the Nextera and Omni buffers were largely interchangeable across all other conditions, while the THS buffer resulted in markedly different profiles in native samples. In-house and commercial enzymes performed similarly. Conclusions We found that the relationship between commonly used measures of library quality differed across temperature and fixation, and so evaluating multiple metrics in assessing the quality of a sample is recommended. Notably, we also found that these choices can bias the functional class of elements profiled and so we recommend evaluating several formulations in any new experiments. Finally, we hope the ATAC-seq workflow formulated in this study on crosslinked samples will help to profile archival clinical specimens.

Recently developed single-cell technologies allow researchers to characterize cell states at ever greater resolution and scale. Caenorhabditis elegans is a particularly tractable system for studying development, and recent single-cell RNA-seq studies characterized the gene expression patterns for nearly every cell type in the embryo and at the second larval stage (L2). Gene expression patterns give insight about gene function and into the biochemical state of different cell types; recent advances in other single-cell genomics technologies can now also characterize the regulatory context of the genome that gives rise to these gene expression levels at a single-cell resolution. To explore the regulatory DNA of individual cell types in C. elegans, we collected single-cell chromatin accessibility data using the sci-ATAC-seq assay in L2 larvae to match the available single-cell RNA-seq data set. By using a novel implementation of the latent Dirichlet allocation algorithm, we identify 37 clusters of cells that correspond to different cell types in the worm, providing new maps of putative cell type–specific gene regulatory sites, with promise for better understanding of cellular differentiation and gene regulation.

We developed an improved assay for single cell profiling of chromatin accessibility that both uses three levels of combinatorial indexing and, in contrast with previous iterations of sci-ATAC-seq and related methods, does not rely on molecularly barcoded Tn5 complexes (sci-ATAC-seq3). Rather, the first two rounds of indexing are achieved by ligation to either end of the conventional, uniformly loaded Tn5 transposase complex (standard Nextera™), while the final round of indexing remains through PCR. Relative to two-level sci-ATAC-seq but similar to sci-RNA-seq3, sci-ATAC-seq3 reduces the per-cell cost of library preparation as well as the rate of collisions, opening the door to experiments on the scale of 10^6 cells. This protocol no longer requires cell sorting, and we also optimized ligase and polymerase choice, kinase concentration, and oligo designs and concentrations, to maximize the number of fragments recovered from each cell. Of note, while maintaining an enrichment in accessible regions, we made the explicit choice to maximize complexity at the expense of specificity for accessible sites. In particular, we found that the fixation conditions could be tuned to adjust the sensitivity (i.e. complexity) vs. specificity (i.e. enrichment in accessible sites) of the assay.

The mechanical microenvironment of primary breast tumors plays a substantial role in promoting tumor progression1. While the transitory response of cancer cells to pathological stiffness in their native microenvironment has been well described2, it is still unclear how mechanical stimuli in the primary tumor influence distant, late-stage metastatic phenotypes across time and space in absentia. Here, we show that primary tumor stiffness promotes stable, non-genetically heritable phenotypes in breast cancer cells. This mechanical memory instructs cancer cells to adopt and maintain increased cytoskeletal dynamics, traction force, and 3D invasion in vitro, in addition to promoting osteolytic bone metastasis in vivo. Furthermore, we established a mechanical conditioning (MeCo) score comprised of mechanically regulated genes as a global gene expression measurement of tumor stiffness response. Clinically, we show that a high MeCo score is strongly associated with bone metastasis in patients. Using a discovery approach, we mechanistically traced mechanical memory in part to ERK-mediated mechanotransductive activation of RUNX2, an osteogenic gene bookmarker and bone metastasis driver3,4. Expanding our conceptual model, we reveal for the first time that a mechanically sensitive transcription factor can maintain a temporal profile of chromatin accessibility in vitro and in vivo. The combination of these distinctive RUNX2 traits permits the stable transactivation of osteolytic target genes that remain upregulated after cancer cells disseminate from their activating microenvironment in order to modify a distant microenvironment. Using genetic, epigenetic, and functional approaches, we could simulate, repress, select and extend RUNX2-mediated mechanical memory and alter cancer cell behavior accordingly. In concert with previous studies detailing the influence of biochemical properties of the primary tumor stroma on distinct metastatic phenotypes5–9, our findings detailing the influence of biomechanical properties support a generalized model of cancer progression in which the integrated properties of the primary tumor microenvironment govern the secondary tumor microenvironment, i.e., soil instructs soil.

The mechanical microenvironment of primary breast tumors plays a substantial role in promoting tumor progression. While the transitory response of cancer cells to pathological stiffness in their native microenvironment has been well described, it is still unclear how mechanical stimuli in the primary tumor influence distant, late-stage metastatic phenotypes across time and space in absentia. Here, we show that primary tumor stiffness promotes stable, non-genetically heritable phenotypes in breast cancer cells. This mechanical memory instructs cancer cells to adopt and maintain increased cytoskeletal dynamics, traction force, and 3D invasion in vitro, in addition to promoting osteolytic bone metastasis in vivo. Furthermore, we established a mechanical conditioning (MeCo) score comprised of mechanically-regulated genes as a global gene expression measurement of tumor stiffness response. Clinically, we show that a high MeCo score is strongly associated with bone metastasis in patients. Using a discovery approach, we mechanistically traced mechanical memory in part to ERK-mediated mechanotransductive activation of RUNX2, an osteogenic gene bookmarker and bone metastasis driver. The combination of these RUNX2 traits permits the stable transactivation of osteolytic target genes that remain upregulated after cancer cells disseminate from their activating microenvironment in order to modify a distant microenvironment. Using genetic, epigenetic, and functional approaches, we were able to simulate, repress, select and extend RUNX2-mediated mechanical memory and alter cancer cell behavior accordingly. In concert with previous studies detailing the influence of biochemical properties of the primary tumor stroma on distinct metastatic phenotypes, our findings detailing the influence of biomechanical properties support a generalized model of cancer progression in which the integrated properties of the primary tumor microenvironment govern the secondary tumor microenvironment, i.e., soil instructs soil.

An improved assay for chromatin accessibility at single-cell resolution in Drosophila melanogaster embryos enables identification of developmental-stage- and cell-lineage-specific patterns of chromatin-level transcriptional regulation. Active gene regulatory elements shape the output of gene transcription and can be mapped across the genome by measuring chromatin accessibility. Eileen Furlong and colleagues apply a technique called ATAC sequencing to profile chromatin accessibility at a single-cell resolution during three stages of Drosophila embryogenesis. They map tissue-specific regulatory elements and show that the chromatin accessibility landscape is sufficient to infer individual cell types and developmental trajectories. A group of cells is found to use regulatory elements of both mesoderm and endoderm, which suggests the existence of a mesendoderm lineage in Drosophila. Understanding how gene regulatory networks control the progressive restriction of cell fates is a long-standing challenge. Recent advances in measuring gene expression in single cells are providing new insights into lineage commitment. However, the regulatory events underlying these changes remain unclear. Here we investigate the dynamics of chromatin regulatory landscapes during embryogenesis at single-cell resolution. Using single-cell combinatorial indexing assay for transposase accessible chromatin with sequencing (sci-ATAC-seq)1, we profiled chromatin accessibility in over 20,000 single nuclei from fixed Drosophila melanogaster embryos spanning three landmark embryonic stages: 2–4 h after egg laying (predominantly stage 5 blastoderm nuclei), when each embryo comprises around 6,000 multipotent cells; 6–8 h after egg laying (predominantly stage 10–11), to capture a midpoint in embryonic development when major lineages in the mesoderm and ectoderm are specified; and 10–12 h after egg laying (predominantly stage 13), when each of the embryo’s more than 20,000 cells are undergoing terminal differentiation. Our results show that there is spatial heterogeneity in the accessibility of the regulatory genome before gastrulation, a feature that aligns with future cell fate, and that nuclei can be temporally ordered along developmental trajectories. During mid-embryogenesis, tissue granularity emerges such that individual cell types can be inferred by their chromatin accessibility while maintaining a signature of their germ layer of origin. Analysis of the data reveals overlapping usage of regulatory elements between cells of the endoderm and non-myogenic mesoderm, suggesting a common developmental program that is reminiscent of the mesendoderm lineage in other species2,3,4. We identify 30,075 distal regulatory elements that exhibit tissue-specific accessibility. We validated the germ-layer specificity of a subset of these predicted enhancers in transgenic embryos, achieving an accuracy of 90%. Overall, our results demonstrate the power of shotgun single-cell profiling of embryos to resolve dynamic changes in the chromatin landscape during development, and to uncover the cis-regulatory programs of metazoan germ layers and cell types.

Objective Type 2 diabetes (T2D) is a complex disease characterized by pancreatic islet dysfunction, insulin resistance, and disruption of blood glucose levels. Genome-wide association studies (GWAS) have identified > 400 independent signals that encode genetic predisposition. More than 90% of associated single-nucleotide polymorphisms (SNPs) localize to non-coding regions and are enriched in chromatin-defined islet enhancer elements, indicating a strong transcriptional regulatory component to disease susceptibility. Pancreatic islets are a mixture of cell types that express distinct hormonal programs, so each cell type may contribute differentially to the underlying regulatory processes that modulate T2D-associated transcriptional circuits. Existing chromatin profiling methods such as ATAC-seq and DNase-seq, applied to islets in bulk, produce aggregate profiles that mask important cellular and regulatory heterogeneity. Methods We present genome-wide single-cell chromatin accessibility profiles in >1,600 cells derived from a human pancreatic islet sample using single-cell combinatorial indexing ATAC-seq (sci-ATAC-seq). We also developed a deep learning model based on U-Net architecture to accurately predict open chromatin peak calls in rare cell populations. Results We show that sci-ATAC-seq profiles allow us to deconvolve alpha, beta, and delta cell populations and identify cell-type-specific regulatory signatures underlying T2D. Particularly, T2D GWAS SNPs are significantly enriched in beta cell-specific and across cell-type shared islet open chromatin, but not in alpha or delta cell-specific open chromatin. We also demonstrate, using less abundant delta cells, that deep learning models can improve signal recovery and feature reconstruction of rarer cell populations. Finally, we use co-accessibility measures to nominate the cell-specific target genes at 104 non-coding T2D GWAS signals. Conclusions Collectively, we identify the islet cell type of action across genetic signals of T2D predisposition and provide higher-resolution mechanistic insights into genetically encoded risk pathways., Highlights • Single cell chromatin accessibility profiles obtained for 1,456 human pancreatic islet cells. • Unique α, β, and δ cell chromatin accessibility signatures identified. • Deep learning approach accurately predicts chromatin accessibility in rare cells. • Nominate effector cell types, causal variants and target genes at diabetes GWAS loci.

ObjectiveType 2 diabetes (T2D) is a complex disease characterized by pancreatic islet dysfunction, insulin resistance, and disruption of blood glucose levels. Genome wide association studies (GWAS) have identified >400 independent signals that encode genetic predisposition. More than 90% of the associated single nucleotide polymorphisms (SNPs) localize to non-coding regions and are enriched in chromatin-defined islet enhancer elements, indicating a strong transcriptional regulatory component to disease susceptibility. Pancreatic islets are a mixture of cell types that express distinct hormonal programs, and so each cell type may contribute differentially to the underlying regulatory processes that modulate T2D-associated transcriptional circuits. Existing chromatin profiling methods such as ATAC-seq and DNase-seq, applied to islets in bulk, produce aggregate profiles that mask important cellular and regulatory heterogeneity.MethodsWe present genome-wide single cell chromatin accessibility profiles in >1,600 cells derived from a human pancreatic islet sample using single-cell-combinatorial-indexing ATAC-seq (sci-ATAC-seq). We also developed a deep learning model based on the U-Net architecture to accurately predict open chromatin peak calls in rare cell populations.ResultsWe show that sci-ATAC-seq profiles allow us to deconvolve alpha, beta, and delta cell populations and identify cell-type-specific regulatory signatures underlying T2D. Particularly, we find that T2D GWAS SNPs are significantly enriched in beta cell-specific and cross cell-type shared islet open chromatin, but not in alpha or delta cell-specific open chromatin. We also demonstrate, using less abundant delta cells, that deep-learning models can improve signal recovery and feature reconstruction of rarer cell populations. Finally, we use co-accessibility measures to nominate the cell-specific target genes at 104 non-coding T2D GWAS signals.ConclusionsCollectively, we identify the islet cell type of action across genetic signals of T2D predisposition and provide higher-resolution mechanistic insights into genetically encoded risk pathways.

SUMMARY While the immediate and transitory response of breast cancer cells to pathological stiffness in their native microenvironment has been well explored, it remains unclear how stiffness-induced phenotypes are maintained over time after cancer cell dissemination in vivo. Here, we show that fibrotic-like matrix stiffness promotes distinct metastatic phenotypes in cancer cells, which are preserved after transition to softer microenvironments, such as bone marrow. Using differential gene expression analysis of stiffness-responsive breast cancer cells, we establish a multigenic score of mechanical conditioning (MeCo) and find that it is associated with bone metastasis in patients with breast cancer. The maintenance of mechanical conditioning is regulated by RUNX2, an osteogenic transcription factor, established driver of bone metastasis, and mitotic bookmarker that preserves chromatin accessibility at target gene loci. Using genetic and functional approaches, we demonstrate that mechanical conditioning maintenance can be simulated, repressed, or extended, with corresponding changes in bone metastatic potential., Graphical Abstract, In brief Watson et al. demonstrate that mechanical conditioning by stiff microenvironments in breast tumors is maintained in cancer cells after dissemination to softer microenvironments, including bone marrow. They show that mechanical conditioning promotes invasion and osteolysis and establish a mechanical conditioning (MeCo) score, associated with bone metastasis in patients.

With the advent of massively parallel sequencing, considerable work has gone into adapting chromosome conformation capture (3C) techniques to study chromosomal architecture at a genome-wide scale. We recently demonstrated that the inactive murine X chromosome adopts a bipartite structure using a novel 3C protocol, termed in situ DNase Hi-C. Like traditional Hi-C protocols, in situ DNase Hi-C requires that chromatin be chemically cross-linked, digested, end-repaired, and proximity-ligated with a biotinylated bridge adaptor. The resulting ligation products are optionally sheared, affinity-purified via streptavidin bead immobilization, and subjected to traditional next-generation library preparation for Illumina paired-end sequencing. Importantly, in situ DNase Hi-C obviates the dependence on a restriction enzyme to digest chromatin, instead relying on the endonuclease DNase I. Libraries generated by in situ DNase Hi-C have a higher effective resolution than traditional Hi-C libraries, which makes them valuable in cases in which high sequencing depth is allowed for, or when hybrid capture technologies are expected to be used. The protocol described here, which involves ∼4 d of bench work, is optimized for the study of mammalian cells, but it can be broadly applicable to any cell or tissue of interest, given experimental parameter optimization.

Genome-wide association studies (GWASs) have become a standard tool for dissecting genetic contributions to disease risk. However, these studies typically require extraordinarily large sample sizes to be adequately powered. Strategies that incorporate functional information alongside genetic associations have proved successful in increasing GWAS power. Following this paradigm, we present the results of 20 different genetic association studies for quantitative traits related to complex diseases, conducted in the Hutterites of South Dakota. To boost the power of these association studies, we collected RNA-sequencing data from lymphoblastoid cell lines for 431 Hutterite individuals. We then used Sherlock, a tool that integrates GWAS and expression quantitative trait locus (eQTL) data, to identify weak GWAS signals that are also supported by eQTL data. Using this approach, we found novel associations with quantitative phenotypes related to cardiovascular disease, including carotid intima-media thickness, left atrial volume index, monocyte count and serum YKL-40 levels.

Single-cell chromatin and RNA analysis Single-cell analyses have begun to provide insight into the differences among and within the individual cells that make up a tissue or organism. However, technological barriers owing to the small amount of material present in each single cell have prevented parallel analyses. Cao et al. present sci-CAR, a pooled barcode method that jointly analyzes both the RNA transcripts and chromatin profiles of single cells. By applying sci-CAR to lung adenocarcinoma cells and mouse kidney tissue, the authors demonstrate precision in assessing expression and genome accessibility at a genome-wide scale. The approach provides an improvement over bulk analysis, which can be confounded by differing cellular subgroups. Science , this issue p. 1380

We applied a combinatorial indexing assay, sci-ATAC-seq, to profile genome-wide chromatin accessibility in ∼100,000 single cells from 13 adult mouse tissues. We identify 85 distinct patterns of chromatin accessibility, most of which can be assigned to cell types, and ∼400,000 differentially accessible elements. We use these data to link regulatory elements to their target genes, to define the transcription factor grammar specifying each cell type, and to discover in vivo correlates of heterogeneity in accessibility within cell types. We develop a technique for mapping single cell gene expression data to single-cell chromatin accessibility data, facilitating the comparison of atlases. By intersecting mouse chromatin accessibility with human genome-wide association summary statistics, we identify cell-type-specific enrichments of the heritability signal for hundreds of complex traits. These data define the in vivo landscape of the regulatory genome for common mammalian cell types at single-cell resolution.

Single cell measurements of gene expression are providing new insights into lineage commitment, yet the regulatory changes underlying individual cell trajectories remain elusive. Here, we profiled chromatin accessibility in over 20,000 single nuclei across multiple stages ofDrosophilaembryogenesis. Our data reveal heterogeneity in the regulatory landscape prior to gastrulation that reflects anatomical position, a feature that aligns with future cell fate. During mid embryogenesis, tissue granularity emerges such that cell types can be inferred by their chromatin accessibility, while maintaining a signature of their germ layer of origin. We identify over 30,000 distal elements with tissue-specific accessibility. Using transgenic embryos, we tested the germ layer specificity of a subset of predicted enhancers, achieving near-perfect accuracy. Overall, these data demonstrate the power of shotgun single cell profiling of embryos to resolve dynamic changes in open chromatin during development, and to uncover thecis-regulatory programs of germ layers and cell types.

Over a million DNA regulatory elements have been cataloged in the human genome, but linking these elements to the genes that they regulate remains challenging. We introduce Cicero, a statistical method that connects regulatory elements to target genes using single cell chromatin accessibility data. We apply Cicero to investigate how thousands of dynamically accessible elements orchestrate gene regulation in differentiating myoblasts. Groups of co-accessible regulatory elements linked by Cicero meet criteria of “chromatin hubs”, in that they are physically proximal, interact with a common set of transcription factors, and undergo coordinated changes in histone marks that are predictive of gene expression. Pseudotemporal analysis revealed a subset of elements bound by MYOD in myoblasts that exhibit early opening, potentially serving as the initial sites of recruitment of chromatin remodeling and histone-modifying enzymes. The methodological framework described here constitutes a powerful new approach for elucidating the architecture, grammar and mechanisms ofcis-regulation on a genome-wide basis.

The extent to which non-coding mutations contribute to Mendelian disease is a major unknown in human genetics. Relatedly, the vast majority of candidate regulatory elements have yet to be functionally validated. Here, we describe a CRISPR-based system that uses pairs of guide RNAs (gRNAs) to program thousands of kilobase-scale deletions that deeply scan across a targeted region in a tiling fashion ("ScanDel"). We applied ScanDel to HPRT1, the housekeeping gene underlying Lesch-Nyhan syndrome, an X-linked recessive disorder. Altogether, we programmed 4,342 overlapping 1 and 2 kb deletions that tiled 206 kb centered on HPRT1 (including 87 kb upstream and 79 kb downstream) with median 27-fold redundancy per base. We functionally assayed programmed deletions in parallel by selecting for loss of HPRT function with 6-thioguanine. As expected, sequencing gRNA pairs before and after selection confirmed that all HPRT1 exons are needed. However, HPRT1 function was robust to deletion of any intergenic or deeply intronic non-coding region, indicating that proximal regulatory sequences are sufficient for HPRT1 expression. Although our screen did identify the disruption of exon-proximal non-coding sequences (e.g., the promoter) as functionally consequential, long-read sequencing revealed that this signal was driven by rare, imprecise deletions that extended into exons. Our results suggest that no singular distal regulatory element is required for HPRT1 expression and that distal mutations are unlikely to contribute substantially to Lesch-Nyhan syndrome burden. Further application of ScanDel could shed light on the role of regulatory mutations in disease at other loci while also facilitating a deeper understanding of endogenous gene regulation.

Conventional methods for profiling the molecular content of biological samples fail to resolve heterogeneity that is present at the level of single cells. In the past few years, single cell RNA sequencing has emerged as a powerful strategy for overcoming this challenge. However, its adoption has been limited by a paucity of methods that are at once simple to implement and cost effective to scale massively. Here, we describe a combinatorial indexing strategy to profile the transcriptomes of large numbers of single cells or single nuclei without requiring the physical isolation of each cell (Single cellCombinatorialIndexing RNA-seq or sci-RNA-seq). We show that sci-RNA-seq can be used to efficiently profile the transcriptomes of tens-of-thousands of single cells per experiment, and demonstrate that we can stratify cell types from these data. Key advantages of sci-RNA-seq over contemporary alternatives such as droplet-based single cell RNA-seq include sublinear cost scaling, a reliance on widely available reagents and equipment, the ability to concurrently process many samples within a single workflow, compatibility with methanol fixation of cells, cell capture based on DNA content rather than cell size, and the flexibility to profile either cells or nuclei. As a demonstration of sci-RNA-seq, we profile the transcriptomes of 42,035 single cells fromC. elegansat the L2 stage, effectively 50-fold “shotgun cellular coverage” of the somatic cell composition of this organism at this stage. We identify 27 distinct cell types, including rare cell types such as the two distal tip cells of the developing gonad, estimate consensus expression profiles and define cell-type specific and selective genes. Given thatC. elegansis the only organism with a fully mapped cellular lineage, these data represent a rich resource for future methods aimed at defining cell types and states. They will advance our understanding of developmental biology, and constitute a major step towards a comprehensive, single-cell molecular atlas of a whole animal.

Sequencing each cell of the nematode Single-cell sequencing is challenging owing to the limited biological material available in an individual cell and the high cost of sequencing across multiple cells. Cao et al. developed a two-step combinatorial barcoding method to profile both single-cell and single-nucleus transcriptomes without requiring physical isolation of each cell. The authors profiled almost 50,000 single cells from an individual Caenorhabditis elegans larval stage and were able to identify and recover information from different, even rare, cell types. Science , this issue p. 661

The extent to which distal non-coding mutations contribute to Mendelian disease remains a major unknown in human genetics. Given that a gene’s in vivo function can be appropriately modeled in vitro, CRISPR/Cas9 genome editing enables the large-scale perturbation of distal non-coding regions to identify functional elements in their native context. However, early attempts at such screens have relied on one individual guide RNA (gRNA) per cell, resulting in sparse mutagenesis with minimal redundancy across regions of interest. To address this, we developed a system that uses pairs of gRNAs to program thousands of kilobase-scale deletions that scan across a targeted region in a tiling fashion (“ScanDel”). As a proof-of-concept, we applied ScanDel to program 4,342 overlapping 1- and 2- kilobase (Kb) deletions that tile a 206 Kb region centered on HPRT1, the gene underlying Lesch-Nyhan syndrome, with median 27-fold redundancy per base. Programmed deletions were functionally assayed by selecting for loss of HPRT1 function with 6-thioguanine. HPRT1 exons served as positive controls, and all were successfully identified as functionally critical by the screen. Remarkably, HPRT1 function appeared robust to deletion of any intergenic or deeply intronic non-coding region across the 206 Kb locus, indicating that proximal regulatory sequences are sufficient for its expression. A sparser mutagenesis screen of the same 206 Kb with individual gRNAs also failed to identify critical distal regulatory elements. Although our screen did find programmed deletions and individual gRNAs with putative functional consequences that targeted exon-proximal non-coding sequences (e.g. the promoter), long-read sequencing revealed that this signal was driven almost entirely by rare, unexpected deletions that extended into exonic sequence. These targeted validation experiments defined a small region surrounding the transcriptional start site as the only non-coding sequence essential to HPRT1 function. Overall, our results suggest that distal regulatory elements are not critical for HPRT1 expression, and underscore the necessity of comprehensive edited-locus genotyping for validating the results of CRISPR screens. The application of ScanDel to additional loci will enable more insight into the extent to which the disruption of distal non-coding elements contributes to Mendelian diseases. In addition, dense, redundant, large-scale deletion scanning with gRNA pairs will facilitate a deeper understanding of endogenous gene regulation in the human genome.

Recent genome-wide association studies (GWAS) have identified a number of novel genetic associations with complex human diseases. In spite of these successes, results from GWAS generally explain only a small proportion of disease heritability, an observation termed the ‘missing heritability problem’. Several sources for the missing heritability have been proposed, including the contribution of many common variants with small individual effect sizes, which cannot be reliably found using the standard GWAS approach. The goal of our study was to explore a complimentary approach, which combines GWAS results with functional data in order to identify novel genetic associations with small effect sizes. To do so, we conducted a GWAS for lymphocyte count, a physiologic quantitative trait associated with asthma, in 462 Hutterites. In parallel, we performed a genome-wide gene expression study in lymphoblastoid cell lines from 96 Hutterites. We found significant support for genetic associations using the GWAS data when we considered variants near the 193 genes whose expression levels across individuals were most correlated with lymphocyte counts. Interestingly, these variants are also enriched with signatures of an association with asthma susceptibility, an observation we were able to replicate. The associated loci include genes previously implicated in asthma susceptibility as well as novel candidate genes enriched for functions related to T cell receptor signaling and adenosine triphosphate synthesis. Our results, therefore, establish a new set of asthma susceptibility candidate genes. More generally, our observations support the notion that many loci of small effects influence variation in lymphocyte count and asthma susceptibility.

Linking regulatory DNA elements to their target genes, which may be located hundreds of kilobases away, remains challenging. Here we introduce Cicero, an algorithm that identifies coaccessible pairs of DNA elements using single-cell chromatin accessibility data, and so connects regulatory elements to their putative target genes. We apply Cicero to investigate how dynamically accessible elements orchestrate gene regulation in differentiating myoblasts. Groups of Cicero-linked regulatory elements meet criteria of "chromatin hubs" — they are enriched for physical proximity, interact with a common set of transcription factors, and undergo coordinated changes in histone marks that are predictive of changes in gene expression. Pseudotemporal analysis revealed that most DNA elements remain in chromatin hubs throughout differentiation. A subset of elements bound by MYOD in myoblasts exhibit early opening in a PBX1- and MEIS1-dependent manner. This strategy can be applied to dissect the architecture, sequence determinants, and mechanisms of cis-regulation on a genome-wide scale.

Plasma lipoprotein(a) (Lp[a]) level is an independent risk factor of cardiovascular disease that is under strong genetic control. We conducted a genome-wide association study of plasma Lp(a) in 386 members of a founder population that adheres to a communal lifestyle, proscribes cigarette smoking, and prepares and eats meals communally. We identified associations with 77 single nucleotide polymorphisms (SNPs) spanning 12.5 Mb on chromosome 6q26-q27 that met criteria for genome-wide significance (Por= 1.3 x 10(-7)) and were within or flanking nine genes, including LPA. We show that variation in at least six genes in addition to LPA are significantly associated with Lp(a) levels independent of each other and of the kringle IV repeat polymorphism in the LPA gene. One novel SNP in intron 37 of the LPA gene was also associated with Lp(a) levels and carotid artery disease number in unrelated Caucasians (P = 7.3 x 10(-12) and 0.024, respectively), also independent of kringle IV number. This study suggests a complex genetic architecture of Lp(a) levels that may involve multiple loci on chromosome 6q26-q27.

Chromatin state and the single cell Identifying the chromatin state of any single cell, which may or may not have a different function or represent different stages relative to others collected within any single culture, experiment, or tissue, has been challenging. Cusanovitch et al. skirted previously identified technological limitations to identify regions of accessible chromatin at single-cell resolution. Combinatorial cellular indexing, a strategy for multiplex barcoding of thousands of single cells per experiment, was successfully used to investigate the genome-wide chromatin accessibility landscape in each of over 15,000 single cells. Science , this issue p. 910

Don't Ape Protein Variation Changes in DNA and messenger RNA (mRNA) expression levels have been used to estimate evolutionary changes between species. However protein expression levels may better reflect selection on divergent and constrained phenotypes. Khan et al. (p. 1100 , published online 17 October; see the Perspective by Vogel ) measured the differences among and within species between mRNA expression and protein levels in humans, chimpanzees, and rhesus macaques, identifying protein transcripts that seem to be under lineage-specific constraint between humans and chimpanzees.

Gene expression differences are shaped by selective pressures and contribute to phenotypic differences between species. We identified 964 copy number differences (CNDs) of conserved sequences across three primate species and examined their potential effects on gene expression profiles. Samples with copy number different genes had significantly different expression than samples with neutral copy number. Genes encoding regulatory molecules differed in copy number and were associated with significant expression differences. Additionally, we identified 127 CNDs that were processed pseudogenes and some of which were expressed. Furthermore, there were copy number-different regulatory regions such as ultraconserved elements and long intergenic noncoding RNAs with the potential to affect expression. We postulate that CNDs of these conserved sequences fine-tune developmental pathways by altering the levels of RNA.

The bacterial composition of the human fecal microbiome is influenced by many lifestyle factors, notably diet. It is less clear, however, what role host genetics plays in dictating the composition of bacteria living in the gut. In this study, we examined the association of ~200K host genotypes with the relative abundance of fecal bacterial taxa in a founder population, the Hutterites, during two seasons (n = 91 summer, n = 93 winter, n = 57 individuals collected in both). These individuals live and eat communally, minimizing variation due to environmental exposures, including diet, which could potentially mask small genetic effects. Using a GWAS approach that takes into account the relatedness between subjects, we identified at least 8 bacterial taxa whose abundances were associated with single nucleotide polymorphisms in the host genome in each season (at genome-wide FDR of 20%). For example, we identified an association between a taxon known to affect obesity (genus Akkermansia) and a variant near PLD1, a gene previously associated with body mass index. Moreover, we replicate a previously reported association from a quantitative trait locus (QTL) mapping study of fecal microbiome abundance in mice (genus Lactococcus, rs3747113, P = 3.13 x 10−7). Finally, based on the significance distribution of the associated microbiome QTLs in our study with respect to chromatin accessibility profiles, we identified tissues in which host genetic variation may be acting to influence bacterial abundance in the gut.

One goal of human genetics is to understand how the information for precise and dynamic gene expression programs is encoded in the genome. The interactions of transcription factors (TFs) with DNA regulatory elements clearly play an important role in determining gene expression outputs, yet the regulatory logic underlying functional transcription factor binding is poorly understood. Many studies have focused on characterizing the genomic locations of TF binding, yet it is unclear to what extent TF binding at any specific locus has functional consequences with respect to gene expression output. To evaluate the context of functional TF binding we knocked down 59 TFs and chromatin modifiers in one HapMap lymphoblastoid cell line. We then identified genes whose expression was affected by the knockdowns. We intersected the gene expression data with transcription factor binding data (based on ChIP-seq and DNase-seq) within 10 kb of the transcription start sites of expressed genes. This combination of data allowed us to infer functional TF binding. Using this approach, we found that only a small subset of genes bound by a factor were differentially expressed following the knockdown of that factor, suggesting that most interactions between TF and chromatin do not result in measurable changes in gene expression levels of putative target genes. We found that functional TF binding is enriched in regulatory elements that harbor a large number of TF binding sites, at sites with predicted higher binding affinity, and at sites that are enriched in genomic regions annotated as “active enhancers.”, Author Summary An important question in genomics is to understand how a class of proteins called “transcription factors” controls the expression level of other genes in the genome in a cell-type-specific manner – a process that is essential to human development. One major approach to this problem is to study where these transcription factors bind in the genome, but this does not tell us about the effect of that binding on gene expression levels and it is generally accepted that much of the binding does not strongly influence gene expression. To address this issue, we artificially reduced the concentration of 59 different transcription factors in the cell and then examined which genes were impacted by the reduced transcription factor level. Our results implicate some attributes that might influence what binding is functional, but they also suggest that a simple model of functional vs. non-functional binding may not suffice.

Epstein-Barr virus (EBV) transformed lymphoblastoid cell lines (LCLs) provide a conveniently accessible and renewable resource for functional genomic studies in humans. The ability to accumulate multidimensional data pertaining to the same individual cell lines, from complete genomic sequences to detailed gene regulatory profiles, further enhances the utility of LCLs as a model system. A lingering concern, however, is that the changes associated with EBV transformation of B cells reduce the usefulness of LCLs as a surrogate model for primary tissues. To evaluate the validity of this concern, we compared global gene expression and methylation profiles between CD20+ primary B cells sampled from six individuals and six independent replicates of transformed LCLs derived from each sample. These data allowed us to obtain a detailed catalog of the genes and pathways whose regulation is affected by EBV transformation. We found that the expression levels and promoter methylation profiles of more than half of the studied genes were affected by the EBV transformation, including enrichments of genes involved in transcription regulation, cell cycle and immune response. However, we show that most of the differences in gene expression levels between LCLs and B cells are of small magnitude, and that LCLs can often recapitulate the naturally occurring gene expression variation in primary B cells. Thus, our observations suggest that inference of the genetic architecture that underlies regulatory variation in LCLs can typically be generalized to primary B cells. In contrast, inference based on functional studies in LCLs may be more limited to the cell lines.