1. Developmental Biology
  2. Genetics and Genomics

Flexible and high-throughput simultaneous profiling of gene expression and chromatin accessibility in single cells

  1. Volker Soltys  Is a corresponding author
  2. Moritz A Peters
  3. Dingwen Su
  4. Marek Kucka
  5. Yingguang Frank Chan  Is a corresponding author
  1. Friedrich Miescher Laboratory of the Max Planck Society, Germany
  2. Department of Translational Genomics, University of Cologne, Germany
  3. University of Groningen, Groningen Institute for Evolutionary Life Sciences, Netherlands
https://doi.org/10.7554/eLife.110034.4
Share this article
Cite this article
  1. Volker Soltys
  2. Moritz A Peters
  3. Dingwen Su
  4. Marek Kucka
  5. Yingguang Frank Chan
(2026)
Flexible and high-throughput simultaneous profiling of gene expression and chromatin accessibility in single cells
eLife 15:RP110034.
https://doi.org/10.7554/eLife.110034.4
  2. Download BibTeX
  3. Download .RIS

eLife Assessment

This study presents a useful methodological advance that better enables the simultaneous measurement of gene expression and chromatin accessibility in individual cells. The evidence supporting the improved detection of gene expression is solid. The method has the potential to be more broadly impactful if it were expanded to include orthogonal validation strategies. This method will be of interest to those studying transcription and gene regulation.

https://doi.org/10.7554/eLife.110034.4.sa0

Abstract

Gene regulation underpins development and is an intricate biological process involving transcription, typically at promoters within accessible chromatin. To understand cell-type-specific regulatory networks, the ability to capture both transcription and chromatin accessibility simultaneously is crucial. However, joint measurements are technically challenging and current methodologies still face adoption challenges. Here, we present easySHARE-seq, an improvement on SHARE-seq for the simultaneous measurement of ATAC- and RNA-seq in single cells. We address several limitations of the previous method by improving the barcode and streamlining the protocol. As a result, easySHARE-seq libraries have a usable sequence of up to 300 bp (+200 bp increase), making it suitable for, e.g., investigation of allele-specific signals or variant discovery. Furthermore, easySHARE-seq libraries do not require a dedicated sequencing run thus saving costs. We applied easySHARE-seq to murine liver nuclei and recovered 19,664 nuclei with joint chromatin and expression profiles. By benchmarking against other combinatorial indexing-based techniques, we showed that we can recover over 1.5-fold more transcripts per cell while retaining high scalability and low cost. To showcase our method, we identified cell types, exploited the multiomic measurements to link cis-regulatory elements to their target genes and investigated liver-specific micro-scale changes. We conclude that easySHARE-seq improves upon previous methods and can produce high-quality multiomic datasets. We expect it to be applicable to a wide range of study designs.

Introduction

Development is an intricate biological process that requires the coordinated expression of thousands of genes, in millions of cells, nearly all of which carry the same genome, at least at the DNA base pair level. The difference between cells therefore emerges through regulation—acting before, during, and after transcription across layers of regulatory mechanism, resulting in lineage specification and proper tissue function (Spitz and Furlong, 2012; Petit et al., 2017). Failures in these mechanisms underlie many human diseases: for instance, activating mutations in the Notch pathway can cause T-cell acute lymphoblastic leukaemia (Weng et al., 2004). There are also many well-documented examples of disrupted enhancer–promoter communication contributing to congenital malformations and cancer (Zabidi and Stark, 2016; Zhang et al., 2013; Spielmann et al., 2012), as well as evolutionary adaptations (Kvon et al., 2016; Senevirathne et al., 2025). Likewise, aberrant DNA methylation reshapes regulatory networks in malignancy (Shen and Laird, 2013). For this reason, functional genomic methods that enable researchers to track developmental processes across regulatory layers are indispensable tools to help us understand gene regulatory processes in health and disease.

Over the past decade, large-scale initiatives such as ENCODE (The ENCODE Project Consortium, 2012) and FANTOM (Forrest et al., 2014) have expanded and spurred the development of a collection of molecular assays to interrogate many regulatory processes, including histone modifications (Allis and Jenuwein, 2016), open chromatin (Buenrostro et al., 2013), (nascent) transcription (Core et al., 2008), and three-dimensional genome conformation (Dekker et al., 2013) (for a detailed review, see Vandereyken et al., 2023). These approaches generally convert functional states into DNA sequence readouts and, together, have transformed our understanding of genome regulation (The ENCODE Project Consortium, 2012). These and other efforts to interrogate and summarise these results have converged on the idea where chromatin accessibility and long-range promoter–enhancer interactions are among the most predictive of gene expression programmes and can bias cell-fate decisions (Petit et al., 2017; Weintraub et al., 2017; Castro et al., 2019; Schoenfelder and Fraser, 2019).

Another inflection point in advances in sequencing and barcoding strategies is the increase in resolution of these functional assays, ultimately being able to distinguish individual cells from each other (‘single-cell’ resolution) (Lähnemann et al., 2020). This change led to the discovery of rare cell populations and dynamic state transitions that were previously masked in conventional bulk assays. Consequently, a range of assays have been adopted to report on chromatin accessibility, transcriptomes, and protein abundance at the single-cell level—namely, scATAC-seq (for chromatin accessibility) (Buenrostro et al., 2015), scRNA-seq (for the transcriptome) (Tang et al., 2009), and scBS-seq (for methylation) (Smallwood et al., 2014). Many of these approaches were based on physical separation between cells, classically in droplets. Paired with a correspondingly large set of DNA barcodes, this made it possible to pool thousands to now millions of cells together in a multiplexed sequencing reaction (Macosko et al., 2015). Yet, because each technology typically captures only one modality, inferring causal relationships between transcriptome and epigenome often relies on computational integration of separately measured datasets, using tools like Seurat (Stuart et al., 2019) or MOFA (Argelaguet et al., 2018). An inevitable assumption across these tools is that there is a shared underlying state across modality. While generally justifiable, it can also obscure fine-grained heterogeneity and complicate attempts at resolving how multiple regulatory layers interact within rare or transient cell types (Clark et al., 2018; Ma et al., 2020).

Thus, the frontier is to measure multiple regulatory layers within the same single cell, enabling direct links between chromatin state and transcriptional output. Several strategies have been developed, each with distinct strengths and drawbacks. Early integrative methods such as scNMT-seq combined chromatin accessibility, DNA methylation, and transcription, but required technically demanding protocols and produced sparse coverage (Clark et al., 2018). Generally, the early dependence on microfluidic chips and instrumentations imposed a limit on usability, throughput, and protocol flexibility. These constraints encouraged the development of alternative approaches. One such approach that offers a simpler workflow and has found broad adoption is combinatorial indexing (Erlich et al., 2009). It relies on successive rounds of a ‘split-and-pool’ procedure such that each cell is likely to be assigned a unique set of barcode segments. The sequential barcode design is powerful for single-cell barcoding as it scales exponentially and the ‘split-and-pool’ procedure is more open to protocol improvements. Earlier examples include sci-CAR, which enabled joint profiling of chromatin accessibility and RNA, yet were limited by modest transcript capture and complex library design (Cao et al., 2018). More recently, paired-seq has been introduced as a paired multiomic approach, increasing capture efficiency for chromatin accessibility while maintaining scalability (Zhu et al., 2019). SHARE-seq improved throughput, robustness, and data quality by refining barcoding strategies for simultaneous chromatin accessibility and transcriptional profiling. However, as with many protocols, its adoption has been constrained by custom sequencing requirements and extremely long barcode lengths (Ma et al., 2020).

Commercial platforms have begun to address some of these barriers. The 10x Genomics Multiome assay offers a streamlined, kit-based solution for parallel scRNA-seq and scATAC-seq with standardised chemistry, but remains expensive and less adaptable to custom study designs. Scale Biosciences and Parse Biosciences provide a fully commercialised combinatorial indexing system that avoids microfluidics but currently does not offer simultaneous capture of scRNA- and scATAC-seq (Brown et al., 2024). Collectively, these methods highlight the central technical challenges of joint profiling: first, balancing sensitivity and breadth across modalities; second, minimising protocol complexity to ensure broad adoption; and last but not least, achieving scalability without compromising data quality. Addressing these limitations remains critical for charting how regulatory programmes interact to specify cell fate and function.

Here, we describe an improved version of SHARE-seq with focus on minimising protocol complexity and maximising applicability. First, we redesigned the barcoding system, focusing on several improvements: (1) Reads can be sequenced up to a length of 300 bp (instead of 100 bp), allowing for variant discovery or assessment of allele-specific signals. (2) All libraries can be multiplexed with standard short-read libraries, decreasing sequencing costs. (3) Decreased protocol length by several hours. (4) The introduction of ‘sub-libraries’ allows for cost-effective sequencing of pilot data. We implemented these changes and applied them to a set of test samples to determine if easySHARE-seq can deliver larger and more flexible experimental designs through its more streamlined protocol and improved sequence quality.

Results

Technical improvements of easySHARE-seq

To develop easySHARE-seq, we implemented several improvements to SHARE-seq (Ma et al., 2020) while retaining the broad overall workflow (Figure 1A). First, we streamlined the barcode design. We retain a combinatorial indexing scheme but shorten each barcode segment to 7 nucleotides (nt), connected by 3 nt overhangs (Hamming distance between barcode segments ≥ 2). This shortens the index sequences of the libraries to 17 nt (from 99 nt in SHARE-seq, Index 1) and 8 nt (Index 2), respectively (Figure 1B, Figure 1—figure supplement 1A; shows a direct comparison). By using two segments of 192 barcodes each plus a separate third segment (to be added in a later PCR step), we achieve a total of 192×192×96 barcode combinations (~3.5 million) in the space of 25 nt. The advantage of this design is that it allows more effective use of sequencing cycles and the sequencing of up to 300 bp of the insert if desired (Figure 1C; compared to 100 nt in SHARE-seq or 90 nt in Chromium Multiome libraries) while keeping the same scalability as SHARE-seq. Having longer insert read lengths increases the recovery of informative DNA variants (Figure 1C), which is crucial in detecting allele-specific signals or cell-specific variant discovery, e.g., in F1 hybrids or cancer cells. Furthermore, this design allows simultaneous sequencing of easySHARE-seq libraries with standard Illumina libraries (‘multiplexing’), which can greatly reduce sequencing costs, especially in smaller pilot experiments. Separately, we shortened and optimised the bench protocol. The single-cell barcoding step is shortened from 4.5 hr to 1.5 hr by replacing all hybridisation steps with a ligation round, as well as eliminating blocking oligos. This has the advantage of minimising library loss while improving yield. Additionally, after barcoding, cells are aliquoted into as many ‘sub-libraries’ as desired, typically 3500 cells each. This results in higher flexibility for the fine-tuning of cell numbers, the cost-effective sequencing of, e.g., a pilot sub-library to assess the overall data quality or the storing of sub-libraries for potential later use. All in all, easySHARE-seq reagents cost USD 0.056/cell, which is similar to the original SHARE-seq protocol. But, we anticipate much lower sequencing costs, as we now circumvent the highly customised 100 nt index reads. For example, even a 50 cycle sequencing kit would result in sufficient information for many count-based applications. A comparison, including cost, between the two protocols, as well as the most commonly used commercial solution, can be found in Table 1. (For a detailed description of the flexibility of easySHARE-seq, instructions on how to modify and incorporate the framework into new designs, how to convert it to a single-channel assay, as well as critical steps to assess when planning to use easySHARE-seq, see Appendix 1.).

Figure 1 with 1 supplement see all
Download asset Open asset
easySHARE-seq enables high-quality and accurate simultaneous scATAC-seq and scRNA-seq profiling.

(A) Schematic workflow of easySHARE-seq. (B) Generation and structure of the cell-specific barcode within Index 1. Total length of the final barcode is just 17 nt compared to 99 nt previously. (C) Fraction of sequenced DNA bases allocated for either barcodes (grey) or genetic information (red) in easySHARE-seq and SHARE-seq using different sequencing kits. (D) Left: Principle of a species-mixing experiment. Murine OP-9 and human HEK cells are mixed prior to easySHARE-seq. After sequencing, sequences associated with each cell barcode are assessed for genome of origin. Middle left: Unique ATAC-seq fragments per cell aligning to the mouse or human genome. Cells are coloured according to their assigned origin (red: human; blue: mouse; orange: doublet). Middle right: Unique RNA-seq transcripts per cell aligning to the mouse or human genome. Right: Percentage of ATAC-seq fragments or RNA-seq transcripts per cell relative to total sequencing reads mapping uniquely to the human genome. 3.17% of all observed cells classified as doublets. Accounting for same-species doublets, this results in a doublet rate of 6.34%. (E) Comparison of unique molecular identifiers (UMIs)/cell across different single-cell technologies. Red shading denotes all multiomic technologies. Datasets are this study, SHARE-seq (Ma et al., 2020) (murine skin cells), sci-CAR (Cao et al., 2018) (murine kidney nuclei), SNARE-seq (Chen et al., 2019) (adult and neonatal mouse cerebral cortex nuclei), 10x Multiome (Bravo González-Blas et al., 2024) (murine liver nuclei), 10x3’ Expression (Su et al., 2021a) (murine liver nuclei), and sci-RNAseq3 (Martin et al., 2023) (E16.5 mouse embryo nuclei). Cells have been downsampled to a common sequencing depth where possible (see Methods). (F) Comparison of fragments per cell across different single-cell technologies. Colouring as in (E). Datasets differing to (A) are 10x scATAC (Nikopoulou et al., 2023) (murine liver nuclei) and sciATAC-seq (Cusanovich et al., 2018) (murine liver nuclei). Cells have been downsampled to a common sequencing depth where possible (see Methods).

Figure 1—source data 1

Unique molecular identifier (UMI) counts per technique for Figure 1E.

https://cdn.elifesciences.org/articles/110034/elife-110034-fig1-data1-v1.zip
Figure 1—source data 2

Fragment counts per technique for Figure 1F.

https://cdn.elifesciences.org/articles/110034/elife-110034-fig1-data2-v1.zip
Figure 1—source data 3

Species mixing counts for Figure 1D.

https://cdn.elifesciences.org/articles/110034/elife-110034-fig1-data3-v1.zip
Figure 1—source data 4

Data for Figure 1C.

https://cdn.elifesciences.org/articles/110034/elife-110034-fig1-data4-v1.xlsx
Table 1
Comparison between multiomic single-cell technologies.
ProtocoleasySHARE-seqSCHARE-seq10x Multiome
Cost/cell0.056$0.043$0.26$
Additional costs?NoneDedicated sequencing runSpecialised instrument
Throughput (cells)>200,000>200,00080,000
Max. sequencing length of insert300 bp100 bp100 bp
Processing time~12 hr~17 hr~14.5 hr
Convertible to single-channel assay?YesNot describedNo

easySHARE-seq labels both transcriptome and accessible chromatin in individual cells

To evaluate the accuracy and cell specificity of the barcoding, we first performed easySHARE-seq on an equally mixed pool between human and murine cell lines (HEK and OP-9, respectively). Such a species-mixing design is a standard way to identify if a given single cell is consistently labelled with a single barcode (Cusanovich et al., 2015). Any barcode that comprises both human and murine reads is indicative of two or more cells sharing the same barcode (‘doublets’; Figure 1D, left). In total, we recovered 3808 cells (see Methods for details on cell filtering). After separately mapping the reads against the human and the mouse transcriptomes and genomes, we examined the proportions of each barcode’s reads that are mapped to the human or the mouse genomes.

Doing this, we found that both chromatin and transcriptome profiles separated well (Figure 1D, middle). Overall, 96.83% of all cells (3687) carried exclusively human or mouse chromatin and transcriptome. Between the modalities, we found that cDNA showed a lower accuracy with increasing transcript counts. This could be the result of lower mapping accuracy among transcripts, as coding regions are more conserved between species compared to intergenic regions (Shabalina and Spiridonov, 2004). Using the above criteria, we identified a total of 124 human–mouse doublets (3.17%; Figure 1D, right). Factoring in undetected human–human and mouse–mouse doublets, we estimated a final doublet rate of 6.34%. For comparison, a 10x Chromium Next GEM experiment with 10,000 cells has a doublet rate of ~7.7% (reported by manufacturer). Importantly, easySHARE-seq doublet rates can easily be lowered further by aliquoting fewer cells within each sub-library, ensuring that there is sufficient barcode diversity to minimise barcode collision. We thus conclude that easySHARE-seq can support single-cell accurately measuring chromatin accessibility and gene expression in single cells.

Performance of easySHARE-seq in murine primary liver cells

To assess the data quality and overall performance of easySHARE-seq, we focused on murine liver. The liver consists of a diverse set of defined primary cell types, ranging from small non-parenchymal cell types such as liver sinusoidal endothelial cells (LSECs) (Aizarani et al., 2019) to large hepatocytes, which potentially may be multinucleated. This provides a balance between tissue diversity and complexity for assessing our protocol. We generated matched high-quality joint chromatin and gene expression profiles for 19,664 adult liver nuclei across four age-matched mice (two male, two female), for an estimated recovery rate of 70.2% (28,000 input cells). On average, we recovered 3629 unique transcripts (Figure 1E, ±2990 transcripts, std) from 1798 genes (±969; Figure 1—figure supplement 1E). We used unique molecular identifiers (UMIs) to label individual transcripts (Figure 1B; see SI for details). For chromatin accessibility, we recovered on average 2213 unique fragments (Figure 1F, ±1995 fragments) from 2011 accessible peaks (±1672 peaks; Figure 1—figure supplement 1E). Closer examination of each data modality showed excellent recovery, which we will discuss in turn. For the transcriptome, 74% of total RNA-seq reads mapped to transcripts in a strand-specific manner (Figure 1—figure supplement 1B). We do note that 69.9% of reads mapped within introns (Figure 1—figure supplement 1C), which is consistent with the use of nuclei for transcriptome profiling. For comparison, this rate is lower than the ~75% intronic reads reported in Ma et al., 2020. For chromatin accessibility, the scATAC-seq fraction of the libraries displayed the characteristic banding pattern during sample preparation (Figure 1—figure supplement 1I), and 55.9% of sequenced fragments are found in peaks on average (Figure 1—figure supplement 1D; ±8.59%, range 30–76.2%). They are also highly enriched at transcription start sites (TSS; mean TSS enrichment score 4.46; Figure 1—figure supplement 1H). For the relationship between sequencing depth and UMI/fragment recovery, see Figure 1—figure supplement 1F, G, and an example track can be seen at Figure 1—figure supplement 1N. For an estimate of ambient RNA contamination, see Figure 1—figure supplement 1O. easySHARE-seq data was also highly reproducible, firstly across sub-libraries (R>0.99) but also between biological replicates (R>0.95; Figure 1—figure supplement 1J, K). These data gave us confidence that easySHARE-seq performed consistently across samples and batches and could form the basis of further in-depth analyses.

We next conducted benchmarking to determine how easySHARE-seq performed relative to other multiomic and representative single-channel assays. To do so, we made use of publicly available datasets, where possible with matched sample type and controlling for sequencing depth by downsampling. Most importantly, we used an identical analytical pipeline to allow proper comparison of results (Figure 1E and F; note that downsampling was not always possible due to lack of raw sequencing data, see Methods for a detailed description and figure legend for tissue type and study). We chose the following multiomic assays for benchmarking: SHARE-seq (Ma et al., 2020), sci-CAR (Cao et al., 2018), SNARE-seq (Chen et al., 2019), and 10x Multiome (Bravo González-Blas et al., 2024). For added context, we also compared against single-modal datasets: for the transcriptome, 10X3’ expression RNA-seq (Su et al., 2021a) and sci-RNAseq3 (Martin et al., 2023); and chromatin, 10x ATACseq (Nikopoulou et al., 2023) and sci-ATACseq (Cusanovich et al., 2018). For the transcriptome, easySHARE-seq recovered significantly higher number of unique transcripts per cell (3629 vs 1183, 1302, 1179, and 2029 compared to SHARE-seq, sci-CAR, SNARE-seq, and 10x Multiome Expression, respectively, p<2×10–16 for all, pairwise t-test, two-tailed). As such, easySHARE-seq performed similarly to single modality assays albeit they still recovered significantly more transcripts per cell (3901 and 4607 for 10x and sciRNAseq3, respectively, p<0.001). Summarised at the level of expressed genes, we also detected a higher number of genes expressed per cell (1488 genes/cell vs 608, 345, 708, and 1222; p<2 × 10–16; Figure 1—figure supplement 1L). Similarly, easySHARE-seq was closer to single-modal assays, even outperforming the 10x3’ expression (1343 for 10x and 2010 for sciRNAseq3, p<0.005). We hypothesise that the higher reported recovery from easySHARE-seq may have benefited from better sample retention as a result of protocol improvements, e.g., having fewer hybridisation-and-wash steps, higher fixation, etc. These improvements could lead to increased UMI counts per gene (better signal dynamic range) and, in turn, more genes passing threshold. However, since we were limited in the number of comparable publicly available datasets, we cannot rule out that the other benchmark datasets may not be fully comparable or representative. Also, despite our best efforts to make a fair comparison, our downsampling procedure may not fully capture the complexity of the other datasets. Nevertheless, seeing that our dataset is ~2 times smaller than the other benchmarking datasets in terms of cells captured, we feel confident to conclude that easySHARE-seq is robust and flexible (via sub-libraries) and may offer an edge over other alternatives. Lastly, it is straightforward to adapt and upscale easySHARE-seq to an scRNA-seq protocol only (see Appendix 1).

Regarding ATAC-seq data quality, easySHARE-seq performed similarly to other published multiomic assays (Figure 1F), but recovered less unique fragments per cell and detected less accessible peaks per cell than SHARE-seq (Figure 1—figure supplement 1M; median fragments in peaks/cell: 1329 [easySHARE-seq] and 3220 [SHARE-seq]). This may be partially due to choosing to prioritise nuclei integrity and RNA-seq quality and therefore performing a higher fixation (0.35% PFA compared to 0.2% in SHARE-seq), which in our experience may result in less efficient transposition reactions. Nevertheless, 94.3% of open chromatin regions identified in this study overlapped those reported in the independent 10x Multiome dataset derived from the same tissue (Bravo González-Blas et al., 2024; Figure 1F), providing cross-dataset validation of our ATAC-seq signal.

Simultaneous scATAC-seq and scRNA-seq profiling in murine primary liver cells

To visualise and identify cell types, we projected the ATAC- and RNA-seq modalities separately into 2D space and subsequently integrated modalities for a combined representation (Hao et al., 2021; Figure 2A, Figure 2—figure supplement 1A, B). Between the modalities, snRNA-seq showed good cluster separation, resulting in eight clusters (Figure 2—figure supplement 1A). Clustering in snATAC-seq alone resulted in only four major clusters, likely due to the lower dynamic range in chromatin accessibility and thus less information content (Figure 2—figure supplement 1B). We then annotated cell types on the integrated dataset based on gene expression of previously established marker genes (Su et al., 2021a; Han et al., 2018). Marker gene expression was highly specific to the clusters (Figure 2B, C, Figure 2—figure supplement 1F) and we recovered all expected cell types (Figure 2—figure supplements 1C, 83% hepatocytes, 7.9% LSECs, 3% hepatic stellate cells [HSCs], 2.04% Kupffer cells, 2.01% B cells, 1.2% neurons, 0.4% monocytes, and 0.3% cholangiocytes). Data recovery per cell differed between cell types, with LSECs having the least UMIs and fragments per cell on average (mean 2644 UMIs/cell and mean 1589 fragments/cell), whereas HSCs had the highest amount of transcripts (mean 3764 UMIs/cell) and neurons the most fragments (mean 3040 fragments/cell; Figure 2—figure supplement 1D and E).

Figure 2 with 1 supplement see all
Download asset Open asset
Cell-type classification in primary liver nuclei using joint expression and chromatin accessibility profiles.

(A) UMAP visualisation of WNN-integrated scRNA-seq and scATAC-seq modalities of 19,664 liver nuclei. Nuclei are coloured by their assigned cell-type identity. (B) WNN-UMAPs of 19,664 liver nuclei with nuclei coloured according to the mean expression strength of marker genes for several cell types (Alb: hepatocytes, Clec4g: liver sinusoidal endothelial cells [LSECs], Dcn: hepatic stellate cells [HSCs], Vsig4: Kupffer cells). Plots for further marker genes can be found in Figure 2—figure supplement 1. (C) Violin plots depicting the distribution of the normalised expression level of marker genes in all cells assigned to each cell type. (D) Pseudo-bulk ATAC-seq coverage tracks of normalised chromatin accessibility in all cells of each cell type in open chromatin regions overlapping marker genes.

Examining our dataset, we also found diversity within cell types. In the liver, hepatocytes fulfil different metabolic functions depending on their physical location in their organising structure called lobules (‘zonation’). These differences in function result in gene expression gradients along zonation, which are well described in hepatocytes (Bravo González-Blas et al., 2024; Paris and Henderson, 2022). This diversity was also reflected in our dataset (Figure 2—figure supplement 1G). For example, zonation markers such as Gls2, Glul, and Cyp3e1 expression showed clear gradients within hepatocytes.

Lastly, examining pseudo-bulk ATAC-seq tracks of the identified cell types showed highly cell-type-specific signals despite lower clustering resolution, indicating that cell-type-specific signals can be identified using each modality independently and showcasing the high congruence between the scATAC-seq and scRNA-seq modalities (Figure 2D). Altogether, our results show that easySHARE-seq generates high quality and reproducible joint cellular profiles of chromatin accessibility and gene expression within primary tissue, expanding our toolkit of multiomic protocols.

Uncovering the cis-regulatory landscape of key regulators through peak–gene associations

The power of multiomic technologies lies in simultaneously capturing information from multiple molecular layers, providing a more holistic view of the underlying biology. For example, connecting a cis-regulatory element (CRE) to its target gene is difficult with separate measurements even though it represents the pivotal step in transcription initiation (Spitz and Furlong, 2012). However, changes to such relationships can impact gene function, and in many cases, cause diseases (Smemo et al., 2014; Lupiáñez et al., 2015; Guenther et al., 2008).

To demonstrate the ability of easySHARE-seq data to support such analyses, we focused on LSECs (1561 total), the cell type with the lowest amount of transcripts and fragments in our dataset (Figure 2—figure supplement 1D, E). Following Ma et al., 2020, we computed the correlation between gene expression and chromatin accessibility at nearby peaks to identify these putative CREs (pCREs, Figure 3A, testing all accessible peaks within ±500 kb of the TSS and controlling for GC content and accessibility strength). In total, we found 81,514 significant peak–gene associations (45% of total peaks, p<0.05) with 15,061 genes having at least one association (65.8% of all genes; Figure 3—figure supplement 1A, B). In rare cases (2.9%), these pCREs were associated with five or more genes. This drops to 0.03% when considering only pCREs within ±50 kb of a TSS (Figure 3—figure supplement 1C, D). These pCREs tended to cluster to regions of higher expressed gene density (2.15 mean expressed genes within 50k bp vs 0.93 for all global peaks), and their associated genes were enriched for biological processes such as mRNA processing, histone modifications, and splicing (Figure 3—figure supplement 1H), possibly reflecting loci with increased regulatory activity. Although potentially informative, these signals may reflect technical artefacts. We therefore restricted the analysis to the most significant gene association per peak, resulting in 40,975 pCREs. Exploring these associations, we found that linked genes exhibited significantly higher expression levels while linked peaks showed greater chromatin accessibility, relative to their non-linked counterparts (Figure 3—figure supplement 1I, J), which could either reflect increased power to detect them or suggest that cis-regulatory associations are enriched at transcriptionally more active loci. Despite this observation, chromatin accessibility and gene expression were uncorrelated within linked pairs (Figure 3—figure supplement 1K; p=0.39, r=0.004). This indicates that chromatin accessibility poorly predicts enhancer activity, consistent with evidence that promoters integrate enhancer signals in a non-linear fashion (Catarino and Stark, 2018; Martinez-Ara et al., 2024).

Figure 3 with 1 supplement see all
Download asset Open asset
Assigning putative cis-regulatory elements (pCREs) to their target genes by correlating simultaneous measurements of gene expression and chromatin accessibility.

(A) Schematic depicting the concept for linking pCREs to their target genes. For each gene, all open chromatin regions (OCRs) within ±500 bp of its transcription start site (TSS) are tested. A pCRE is linked to a gene if the Spearman correlation of its chromatin accessibility to the gene’s expression falls outside the expected distribution estimated by correlating chromatin accessibility of 100 unrelated peaks (on different chromosomes) to the gene expression. (B) Genes ranked by their number of significantly correlated pCREs (p<0.05, ±500 kbp from TSS) in liver sinusoidal endothelial cells (LSECs). Marked are genes in the top 1% that are either transcription factors or cell-type-specific regulators shown to fulfil a critical role in LSECs. (C) pCREs are enriched for TSS proximity. Normalised density of all open chromatin regions within ±50 kbp of a TSS (red) and of all pCREs within ±50 kbp of a TSS (blue). (D) Aggregate snATAC-seq pileup (red) of LSECs at the Gata4 locus and 500 kbp upstream region. Grey bars indicate open chromatin regions. Loops denote pCREs significantly correlated with Gata4 and are coloured by Spearman correlation of respective pCRE–Gata4 comparison.

Figure 3—source data 1

Peak–gene associations.

https://cdn.elifesciences.org/articles/110034/elife-110034-fig3-data1-v1.zip

We then ranked genes based on their number of associated pCREs (Figure 3B). Within the top 1% of genes with the most pCRE associations were many key regulators and transcription factors. Examples include Taf5, which directly binds the TATA-box (Chen et al., 2021) and is required for the initiation of transcription, or Gata4, which has been identified as the master regulator for LSEC specification during development. For instance, LSEC-specific knock-out of Gata4 in adult or embryonic mice both lead to transdifferentiation of discontinuous LSEC into continuous capillaries, resulting in liver hypoplasia and fibrosis, even lethality in the case of embryonic conditional loss-of-function (Winkler et al., 2021; Géraud et al., 2017). In adult mice, LSECs also control regeneration and metabolic maturation of liver tissue in adult mice. As such, it incorporates a variety of signals and its expression needs to be strictly regulated, which can be related to its numerous pCRE associations (eight total; Figure 3D). Similarly, Igf1 also integrates signals from many different pCREs (Lara-Diaz et al., 2017) (nine total; Figure 3—figure supplement 1G). Notably, pCREs are significantly enriched at TSS, even relative to background enrichment in global peaks (Figure 3C). We thus conclude that easySHARE-seq makes it possible to directly investigate the relationship between chromatin accessibility and gene expression and link pCREs to their target genes at genomic scale, even in relatively rare cell types with low mRNA contents.

De novo identification of open chromatin regions and genes displaying zonation in LSECs

Micro-scale changes in gene regulation or expression can often have a major impact, not only within a cell through cell-fate determination but also via cell–cell communications, for example during embryonic development or in certain diseases (Kvon et al., 2016; Horn et al., 2013). To demonstrate the ability of easySHARE-seq to capture and assess these micro-scale changes, we investigated the process of zonation in LSECs. The liver consists of hexagonal units called lobules where blood flows from the portal vein and arteries towards a central vein (Baratta et al., 2009; Jungermann and Kietzmann, 1996; Figure 4A). The central–portal (CP) axis is characterised by a morphogen gradient, e.g., Wnt2, secreted by central vein LSECs and thus activating the canonical Wnt pathway, giving rise to spatial division of labour among cells along it (‘zonation’) (Braeuning et al., 2006; Planas-Paz et al., 2016; Wang et al., 2015). While studying zonation in hepatocytes is more common, doing so in non-parenchymal cells such as LSECs is challenging as these are physically small cells with low mRNA content (Figure 2—figure supplement 1D, E), often lying below the detection limit of current spatial transcriptomic techniques. As a result, only very few studies assess zonation in LSECs on a genomic level (Halpern et al., 2018). However, LSECs are critical to liver function as they line the artery walls, clear, and process endotoxins. As outlined above in discussing its master regulator Gata4, LSECs are necessary for proper liver regeneration and setting up morphogen gradients to regulate hepatocyte gene expression (Knolle and Wohlleber, 2016; Smedsrød, 2004; Rafii et al., 2016). For these reasons, a fine-grained understanding of these cells and their gene regulation is crucial for tackling many diseases.

Figure 4
Download asset Open asset
Zonation profiles in liver sinusoidal endothelial cells (LSECs) across gene expression and chromatin accessibility.

(A) Schematic depiction of a liver lobule. A liver lobule has a ‘central–portal axis’ defined by morphogen gradients starting from the central vein to the portal vein and portal artery. The sinusoidal capillary channels are lined with LSECs. (B) UMAP of 1561 LSECs coloured by pseudotime. (C) Changes along the central–portal axis at the Wnt2 locus. Top: Aggregate snATAC-seq profile (red) of LSECs at the Wnt2 locus. Grey bars denote identified open chromatin regions (OCRs). Bottom: In blue, LOESS trend line of mean normalised Wnt2 gene expression along pseudotime/the central–portal axis (CP axis, central vein, CV; portal vein, PV). In red, LOESS trend line of mean normalised chromatin accessibility in OCRs at the Wnt2 locus along the CP axis. (D) LOESS trend line of mean normalised gene expression (blue) of marker genes and mean normalised chromatin accessibility (red) at OCRs overlapping the marker gene along the CP axis for pericentral markers (top, increased towards the central vein, Kit, Dkk3, and Rspo3) and periportal markers (increased towards the portal vein, Efnb2, Meis1, and Ltbp4). (E) Left: Zonation profiles of 550 genes along the CP axis. Right: Zonation profiles of 744 open chromatin regions along the CP axis. All profiles are normalised by their maximum.

Figure 4—source data 1

Marker genes across zonation for Figure 4E.

https://cdn.elifesciences.org/articles/110034/elife-110034-fig4-data1-v1.zip
Figure 4—source data 2

Open chromatin regions across zonation for Figure 4E.

https://cdn.elifesciences.org/articles/110034/elife-110034-fig4-data2-v1.zip

We therefore asked if we can recover zonation gradients and potentially identify novel marker genes and open chromatin regions displaying zonation from our data. Specifically, we ordered LSECs along pseudotime and then checked marker gene expression and chromatin accessibility gradients along it, asking if this recapitulates zonation gradients discovered using paired-cell sequencing (Figure 4B–D; Halpern et al., 2018). This analysis recovered the expected zonation profiles for several known zonation markers. For example, Wnt2 expression decreased strongly along the CP axis as did chromatin accessibility of all three peaks at the Wnt2 locus (Figure 4C). We also recovered the expected zonation profiles for several further known pericentral (decrease along the CP axis, e.g. Rspo3, another canonical Wnt agonist and Kit) and periportal (increase along the CP axis, e.g. Efnb2 and Meis1) marker genes, as well as chromatin accessibility gradients at their associated open chromatin regions (Figure 4D). Knowing that the imputed pseudotime faithfully ordered cells along their zonation gradients, we sought to identify novel genes and open chromatin regions displaying zonation in LSECs based on the decrease or increase of the rolling mean in expression or accessibility along pseudotime (see Methods). In total, we classified 179 genes and 198 open chromatin regions as pericentral, and 371 genes and 546 open chromatin regions showed periportal zonation profiles (Figure 4E). The list of markers contained many genes regulating epithelial growth and angiogenesis (e.g. Efna1, Nrg2, Jag1, Bmp2, Bmper, Nrp2) (Theilmann et al., 2020; Russell et al., 1999; Vihanto et al., 2005; Duan et al., 2018), related to regulating hepatocyte functions and communication (e.g. Dpp4, Foxo1, Foxp1, Insr, Ephb4) (Dong et al., 2008; Das et al., 2010; Ghorpade et al., 2018; Zou et al., 2015), as well as immunological functions (e.g. Il1a, Pigr, Tgfb1) (Dewidar et al., 2019; Barbier et al., 2014), suggesting that these processes show variation along the PC axis. As many of these newly identified markers are implicated in liver illnesses such as cirrhosis, fibrosis, or non-alcoholic fatty liver disease (Miyao et al., 2015; Su et al., 2021b), these genes are potential new biomarkers for their identification and the open chromatin regions starting points for investigating the role of gene regulation in their emergence. This demonstrates the sensitivity and applicability of easySHARE-seq to detect micro-scale changes even in lowly expressed cell types and simultaneously underscores the advantage of multiomic measurements, where pseudotime inferred from gene expression can then be linked to chromatin state dynamics.

Discussion

Understanding complex processes such as gene regulation or disease states requires the integration of multiple layers of information. To this end we have developed easySHARE-seq to streamline the generation of high-quality joint profiling of chromatin accessibility and gene expression within single cells in a way that can accommodate diverse experimental designs. In this section, we will discuss the following topics in turn: the improvements, the implications on experimental designs, current limitations, and further improvements.

Protocol improvements

We noted that despite the substantial impact from the original SHARE-seq publication, there has not been widespread adoption of the protocol itself by other groups. One possible explanation could be the sheer number of customisations required (ordering and building the original barcodes, custom sequencing runs, etc.), which poses a high barrier of entry for many other users. Taking these into account, we have incorporated a number of improvements in developing easySHARE-seq: barcode simplification, reallocation of read lengths, protocol streamlining, and the ability to split a larger sample into ‘sub-libraries’.

The most consequential change in the protocol was our redesign of the barcodes. Here, we have optimised for barcode complexity within a short sequence space, while preserving sufficient molecular stability throughout the multiple rounds of handling and buffer changes. The species-mixing experiment, together with the high consistency in cell-type-specific signals, gave us confidence that our strategy was sufficiently robust for multi-modal single-cell profiling. Having significantly shortened the barcodes (from 107 nt to 25 nt), we could now retain much of the read lengths for the transcripts or fragments themselves (in other words: data). We argue that just read counts in transcriptome or chromatin accessibility profiling capture only some—but far from all—aspects of gene regulation. There are entire aspects of regulatory mechanisms that benefit from longer read lengths such as variant detection and subsequent analysis of allele-specific signals, RNA-editing, insertion/deletions in CREs, and more. These applications may be entirely missed under the original SHARE-seq configuration. Next, a number of our protocol improvements, when taken in combination, go beyond their individual effects. For instance, by shortening the protocol (12 vs 17 hr) and streamlining various enzymatic steps, we showed that we can achieve better sensitivity in RNA-seq (greater UMIs per cell, among others; Figure 1E). This, together with our introduction of sub-libraries using the i5 barcoding segment, enabled both higher throughput and greater flexibility in experimental design. This, in turn, allows users to conduct smaller-scale pilot tests, while retaining the option to increase sequencing efforts after confirming the quality of the libraries. In our experience, such flexibility is crucial in day-to-day experimental designs and is often not captured or not possible in previous protocols. In terms of costs per cell, easySHARE-seq performs similarly to standard SHARE-seq with ~5.6 cents/cell, which is before factoring in the reduction in sequencing costs compared to standard SHARE-seq, and a fraction of the costs (<25%) of commercially available platforms without considering their specialised instrument costs.

Experimental designs

As alluded to in the last section, we think easySHARE-seq holds a number of technical advantages that translate well into an ability to support more (and different) experimental designs. For example, having the longer read-lengths and thus better power to resolve genomic variants will help in the majority of cases where the samples are diverse, from non-inbred individuals or may carry de novo mutations as in cancer. Having sub-libraries in combination with shorter experimental time makes easySHARE-seq more practical when performing multiple experiments, e.g., assaying more conditions or biological replicates, as in practice, one big experiment with hundreds of thousands of cells rarely reflects actual experimental design.

Our results from the murine LSECs also show that we can recover cells that may be rare or have very low mRNA content. This, in turn, implies that we may be able to use easySHARE-seq when investigating diseases originating from rare cell populations. For example, pancreatic β-cells make up only 1–2% of the endocrine tissue in the pancreas, yet their failure to function has major implications for human health (Dludla et al., 2023). As demonstrated in the example of LSEC zonation and Gata4, the key cell–cell signalling hub during cell specification may involve specific cell populations to set up morphogen gradients like Wnt2, or specific genes like Gata4 to initiate regulatory cascades. In cases like this, researchers may benefit from having access to both regulatory channels, as solely focusing on gene expression or chromatin accessibility alone may not reveal the true significance of the shift. Another area where this is potentially beneficial is in evolutionary studies. For example, this allows exploring how gene-enhancer dynamics change between species or populations and how that in turn might shape differences in the transcriptome.

Limitations and future developments

While we are encouraged by the improvements we have implemented in easySHARE-seq, there are still a number of limitations that will require future work to address. For example, single-channel assays still produce higher quality data compared to our multiomic protocol. Additionally, validating this protocol in a widely used cell line would provide independent confirmation of its improvements. Second, there is opportunity to increase ATAC-seq data quality further without compromising the RNA-seq channel. Compared with SHARE-seq, however, it somewhat lacks resolution. Potential experimental steps that can be optimised are different aspects of fixation or tagmentation; for a detailed description, see Appendix 1. We acknowledge that the current relatively lower quality of the ATAC-seq data may introduce increased variance in downstream analyses, particularly in the context of multi-modal integration with gene expression data making data interpretation more challenging. Lastly, while easySHARE-seq is significantly less expensive than 10x or SHARE-seq (factoring in sequencing costs), it still requires a significant initial investment for the DNA oligos.

In terms of future developments of this assay, other potential improvements might be the introduction of sample-specific barcodes to allow multiplexing many samples at once. While in the RNA-seq, this can be achieved rather easily (via the RT primer, see Appendix 1), this is more challenging in the ATAC-seq. Furthermore, there is potential to decrease the amount of reagents, such as T4 ligase, thus decreasing overall costs and increasing barcode complexity further to handle millions of cells simultaneously.

Conclusion

In conclusion, we present here easySHARE-seq, including benchmarking against existing methods and real-life examples. While the presented biological examples may be more limited in scope, it contains a number of cell types, including examples of micro-differentiation that goes some way towards demonstrating that easySHARE-seq is not only applicable to a variety of study designs but also that the simultaneous measurements can be used to disentangle the complex molecular hierarchy and relationships between different regulatory layers. We envision easySHARE-seq as an important technological step towards the more widespread and diverse use of multiomic techniques. Ultimately, these techniques will play an important role in understanding gene regulation in health and disease, differentiation or lineage commitment, and determining genetic variants affecting those processes.

Methods

Key resources table
Reagent type (species) or resourceDesignationSource or referenceIdentifiersAdditional information
Strain, strain background (Mus musculus)C57BL/6JCharles River LaboratoriesRRID:IMSR_JAX:000664
Strain, strain background (M. musculus)PWD/PhJJackson LaboratoriesRRID:IMSR_JAX:004660
Cell line (Homo sapiens)Embryonic kidney cells HEK293RRID:CVCL_0045
Cell line (M. musculus)OP9 Stromal cellsGift from Juan Carlos Zúñiga-PflückerRRID:CVCL_4398
Sequence-based reagentTn5-APicelli et al., 2014Transposase oligoTCGTCGGCAGCGTCAGATGTGTATAAGAGACAG
Sequence-based reagentTn5-B2SPicelli et al., 2014Transposase oligoCAGGGCATTTGCACCCCATGCC
Sequence-based reagentTn5-reversePicelli et al., 2014Transposase oligo/5Phos/CTGTCTCTTATACACATCT
Sequence-based reagentReverse-Transcription PrimerThis paperReverse-Transcription primer/5Phos/GGGCTCGGAGATGTGTATAAGAGACAGNNNNNNNNNNT/iBiodT/TTTTTTTTTTTTTTTTTTTTTTTTVN
Sequence-based reagentTS oligoThis paperTemplate Switch OligorUrUrUTCGTCGGCAGCGTCAGATGTGTATAAGAGACArGrGrG
Sequence-based reagentI7-Tru-Seq-long primerThis paperI7-Tru-Seq-long primerCAAGCAGAAGACGGCATACGAGAT
Recombinant ProteinMaxima-H Reverse TranscriptaseThermo FisherEP0751cDNA generation, template switching
Recombinant ProteinTn5This paperATAC-seq, library generation
Recombinant ProteinSUPERaseThermo FisherAM2694RNase Inhibitor
Recombinant ProteinRecombinant RNase InhibitorJena BiosciencePCR-392SRNase Inhibitor
Recombinant ProteinT4 LigaseNew England BiolabsM0202LLigation during barcoding
Recombinant ProteinQ5 PolymeraseNew England BiolabsM0491LPCR amplification
Recombinant ProteinProteinase KITWA4392,0010Reverse crosslinking
Recombinant ProteinDynabeads M280Thermo Fisher112.05DcDNA pull-down
Recombinant ProteinRNaseInQIAGENY9240LTemplate switching
Chemical compound, drugDAPIThermo FisherD1306Cell counting
Chemical compound, drugFormaldehydeThermo Fisher28906Fixation
Chemical compound, drugNP-40Sigma-Aldrich492016Permeabilisation
Chemical compound, drugTween-20Bio-Rad#1706531Permeabilisation
Chemical compound, drugTAPSCarl Roth6982.4Buffer
Chemical compound, drugDMFCarl Roth6251.1Buffer
Chemical compound, drugPEG6000Sigma-Aldrich528877Buffer
Chemical compound, drugSize Selection BeadsSpeed BeadsSelection of DNA fragments
Commercial assay or kitMinElute PCR Cleanup KitQIAGEN28004Purification
Commercial assay or kitQubitThermo FisherQ33230Quantification of DNA
SoftwarebwaLi and Durbin, 2009bwa v.0.7.17, RRID:SCR_010910Read alignment
SoftwareUMI-toolsSmith et al., 2017UMI-tools v.1.1.2, RRID:SCR_017048Processing RNA-seq
SoftwarefeatureCountsLiao et al., 2014featureCounts v.2.0.1, RRID:SCR_012919Counting transcripts
SoftwareSTARDobin et al., 2013
STAR v.2.7.9a, RRID:SCR_004463Aligning RNA-seq
SoftwareSeuratHafemeister and Satija, 2019Seurat v. 5.0.1, RRID:SCR_016341Analysis RNA-seq
SoftwareSignacStuart et al., 2021Signac v.1.12.9, RRID:SCR_021158Analysis ATAC-seq
SoftwareMACS2Zhang et al., 2008macs2 v.2.2.7.1, RRID:SCR_013291Peak calling
SoftwareMonocleTrapnell et al., 2014monocle v.2.22, RRID:SCR_018685Pseudotime
SoftwareclusterProfiler – GO enrichmentYu et al., 2012clusterProfiler v3.18.1, RRID:SCR_016884GO enrichment

Animal model and tissue preparation

Mice

All animal experimental procedures were carried out under the licence number EB 01-21M at the Friedrich Miescher Laboratory of the Max Planck Society in Tübingen, Germany. The procedures were reviewed and approved by the Regierungspräsidium Tübingen, Germany. Liver was collected from both male and female wild-type C57BL/6 and PWD/PhJ mice aged between 9 and 11 weeks.

Study design

From each strain, we generated easySHARE-seq libraries for one male and one female mouse (four total). For each individual, we sequenced two sub-libraries, resulting in eight easySHARE-seq libraries.

Cell culture

For the species-mixing experiment, HEK cells were cultured in media containing DMEM/F-12 with GlutaMAX Supplement, 10% FBS, and 1% Penicillin-Streptomycin (PenStrep) at 37°C and 5% CO2. Cells were harvested on the day of the experiment by simply pipetting them off the plate and were then spun down for 5 min at 250G.

Murine OP9-DL4 cells were a gift from Juan Carlos Zúñiga-Pflücker. They were cultured in alpha-MEM medium containing 5% FBS and 1% PenStrep. On the day of the experiment, the cells were harvested by aspirating the medium and adding 4 ml of trypsin, followed by an incubation at 37°C for 5 min. Then, 5 ml of medium was added and cells were spun down for 5 min at 250G. After counting both cell lines using TrypanBlue and the Evos Countess II, equal cell numbers were mixed.

Cell lines were periodically tested for mycoplasma status.

Liver nuclei

The liver was extracted, rinsed in HBSS, cut into small pieces, frozen in liquid nitrogen, and stored in the freezer at –80°C for a maximum of 2 weeks. On the day of the experiment, 1 ml of ice-cold Lysis Solution (0.1% Triton X-100, 1 mM DTT, 10 mM Tris-HCl pH 8, 0.1 mM EDTA, 3 mM Mg(Ac)2, 3 mM CaCl2, and 0.32 M sucrose) was added to the tube. The cell suspension was transferred to a pre-cooled Douncer and dounced 10x using Pestle A (loose) and 15x using Pestle B (tight). The solution was added to a thick wall ultracentrifuge tube on ice and topped up with 4 ml ice-cold Lysis Solution. Then 9 ml of sucrose solution (10 mM Tris-HCl pH 8.0, 3 mM Mg(Ac)2, 3 mM DTT, 1.8 M sucrose) was carefully pipetted to the bottom of the tube to create a sucrose cushion. Samples were spun in a pre-cooled ultracentrifuge with an SW-28 rotor at 24,400 rpm for 1.5 hr at 4°C. Afterwards, all supernatant was carefully aspirated so as not to dislodge the pellet at the bottom and 1 ml ice-cold DEPC-treated water supplemented with 10 µl SUPERase and 15 µl Recombinant RNase Inhibitor was added. Without resuspending, the tube was kept on ice for 20 min. The pellet was then resuspended by pipetting ~15 times slowly up and down followed by a 40 µm cell straining step. Counting of the nuclei using DAPI and the Evos Countess II was immediately followed up by fixation.

easySHARE-seq protocol

Preparing the barcoding oligonucleotides

There are two barcoding rounds in easySHARE-seq with 192 unique barcodes distributed across two 96-well plates in each round (see Supplementary file 1 for a full list of oligonucleotide sequences). Each barcode (BC) is pre-annealed as a DNA duplex for improved stability. The first round of barcodes contains two single-stranded linker sequences at its ends, as well as a 5’ phosphate group to ligate the different barcodes together. The first single-stranded overhang links the barcode to a complementary overhang at the 5’ end of the cDNA molecule or transposed DNA molecule, which originates from either the RT primer or the Tn5 adapter. The second overhang (3 bp) is used to ligate it to the second round of barcodes (Figure 1B). Each duplex needs to be annealed prior to cellular barcoding, preferably on the day of the experiment. No blocking oligos are needed.

The Round1 BC plates contain 10 µl of 4 µM duplexes in each well and Round2 BC plates contain 10 µl of 6 µM barcode duplexes in each well, all in Annealing Buffer (10 mM Tris pH 8.0, 1 mM EDTA, 30 mM KCl). Pre-aliquoted barcoding plates can be stored at –20°C for at least 3 months. On the day of the experiment, the oligo plates were thawed and annealed by heating plates to 95°C for 2 min, followed by cooling down the plates to 20°C at a rate of –2°C per minute. Finally, the plates were spun down. Until the annealed barcoding plates are needed, they should be kept on ice or in the fridge.

This barcoding scheme is very flexible and currently supports a throughput of ~350,000 cells (assuming 96 indexing primers) per experiment, limited only by sequencing cost and availability of indexing primer. The barcodes were designed to have at least a Hamming distance of 2. See Appendix 1 for further details on the barcoding system and flexibility.

Tn5 preparation

Tn5 was expressed in-house as previously described (Picelli et al., 2014). Two differently loaded Tn5 are needed for easySHARE-seq, one for the tagmentation, loaded with an adapter for attaching the first barcodes (termed Tn5-B2S), and one for library preparation with a standard Illumina sequencing adapter (termed Tn5-A-only). See Supplementary file 1 for all sequences.

To assemble Tn5-B2S, two DNA duplexes were annealed: 20 µM Tn5-A oligo with 22 µM Tn5-reverse and 20 µM Tn5-B2S with 22 µM Tn5-reverse, all in 50 mM NaCl and 10 mM Tris pH 8.0. Oligos were annealed by heating the solution to 95°C for 30 s and cooling it down to 20°C at a rate of 2 °C/min. An equal volume of duplexes was pooled and then 200 µl of unassembled Tn5 was mixed with 16.5 µl of duplex mix. The Tn5 was then incubated at 37°C for 1 hr, followed by 4°C overnight. The Tn5 can then be stored at –20°C. In our hands, Tn5 did not show a decrease in activity after 10 months of storage.

To assemble Tn5-A-only, 10 µM of Tn5-A and 10.5 µM Tn5-reverse was annealed using the same conditions as described above. Again, 200 µl of unassembled Tn5 was mixed with 16.5 µl of Tn5-A duplex and incubated at 37°C for 1 hr, followed by 4°C overnight. The Tn5 can then be stored for later and repeated use for more than 10 months at –20°C. We observed an increase in all Tn5 activity during the first months of storage, possibly due to continued transposome assembly in storage.

Fixation

One million liver nuclei (‘cells’ for short) were added to ice-cold PBS for 4 ml total. After mixing, 87 µl 16% formaldehyde solution (0.35%; for liver nuclei) or 25 µl 16% formaldehyde solution (0.1%; for HEK and OP9 cells) was added and the suspension was mixed by pipetting up and down exactly three times with a P1000 pipette set to 700 µl. The suspension was incubated at room temperature for 10 min. Fixation was stopped by adding ice-cold Stop-Mix (224 µl 2.5 M glycine, 200 µl 1 M Tris-HCl pH 8.0, 53 µl 7.5% BSA in PBS). The suspension was mixed exactly three times with a P1000 pipette set to 850 µl and incubated on ice for 3 min followed by centrifugation at 500G for 5 min at 4°C. Supernatant was removed and the pellet was resuspended in 1 ml Nuclei Isolation Buffer (NIB; 10 mM Tris pH 8.0, 10 mM NaCl, 2 mM MgCl2, 0.1% NP-40) and kept on ice for 3 min followed by straining the suspension with a 40 µm cell strainer. It was then spun down at 500G for 3 min at 4°C and resuspended in ~100–200 µl PBSi (1x PBS+0.4 U/µl Recombinant RNaseInhibitor, 0.04% BSA, 0.2 U/µl SUPERase, freshly added), depending on the amount of input cells. Cells were then counted using DAPI and the Countess II and concentration was adjusted to 2 M cells/ml using PBSi.

Tagmentation

In a typical easySHARE-seq experiment for this study, eight tagmentation reactions with 10,000 cells each followed by three total RT reactions were performed. This results in sequencing libraries for around 30000 cells. To increase throughput, simply increase the amount of tagmentation and RT reactions accordingly. No adjustment is needed to the barcoding. Each tube and PCR strip until the step of Reverse Crosslinking was coated before use by rinsing it with PBS+0.5% BSA to maximise cell recovery. Even though all cells could theoretically be tagmented or reverse-transcribed in one big reaction, we found separating them into smaller volumes resulted in better data quality.

For each tagmentation reaction, 5 µl of 5x TAPS buffer, 0.25 µl 10% Tween, 0.25 µl 1% Digitonin, 3 µl PBS, 1 µl Recombinant RNaseInhibitor, and 9 µl of H2O were mixed. TAPS buffer was made by first making a 1 M TAPS stock solution in H2O, followed by adjustment of the pH to 8.5 by titrating 10 M NaOH. Then, 4.25 ml H2O, 500 µl 1 M TAPS pH 8.5, 250 µl 1 M MgCl2, and 5 ml N-N-dimethylformamide (DMF) were mixed on ice and in order. When adding DMF, the buffer heats up so it is important to be kept on ice. The resulting 5x TAPS buffer can then be stored at 4°C for short-term use (1–2 months) or for long-term storage at –20°C (>6 months). Then, 5 µl of cell suspension at 2 M cells/ml in PBSi was added to the tagmentation mix for each reaction, mixed thoroughly, and finally 1.5 µl of Tn5-B2S was added. The reaction was incubated on a shaker at 37°C for 30 min at 850 rpm. Afterwards, all reactions were pooled on ice into a pre-cooled 15 ml tube. The reaction wells were washed with ~30 µl PBSi, which was then added to the pooled suspension in order to maximise cell recovery. The suspension was then spun down at 500G for 3 min at 4°C. Supernatant was aspirated and the cells were washed with 200 µl NIB followed by another centrifugation at 500G for 3 min at 4°C.

We only observed cell pellets when centrifuging after fixation and only when using cell lines as input material. Therefore, when aspirating supernatant at any step, it is preferable to leave around 20–30 µl liquid in the tube. Additionally, it is recommended to pipette gently at any step as not to damage and fracture the cells.

Reverse transcription

As stated above, three tagmentation reactions were combined into one RT reaction. When increasing cells to more than 30,000 per RT reaction, we observed a steep drop in reaction efficiency. The Master Mix for one RT reaction contained 3 µl 100 µM RT primer, 2 µl 10 mM dNTPs, 6 µl 5x MaximaH RT Buffer, 4.5 µl 50% PEG6000, 1.5 µl H2O, 1.5 µl SUPERase, and 1.66 µl MaximaH RT. The RT primer contains a polyT tail, a 10 bp UMI sequence, a biotin molecule, and an adapter sequence used for ligating onto the first round of barcoding oligos.

The cell suspension was resuspended in 10 µl NIB per RT reaction and added to the Master Mix for a total of 30 µl. As PEG is present, it is necessary to pipette ~30 times up and down to ensure proper mixing. The RT reaction was performed in a PCR cycler with the following protocol: 52°C for 12 min; then 2 cycles of 8°C for 12 s, 15°C for 45 s, 20°C for 45 s, 30°C for 30 s, 42°C for 2 min, and 50°C for 3 min. Finally, the reaction was incubated at 52°C for 5 more minutes. All reactions were then pooled on ice into a pre-cooled and coated 15 ml tube and the reaction wells were washed with ~40 µl NIB, which was then added to the pooled cell suspension in order to maximise cell recovery. The suspension was then spun down at 500G for 3 min at 4°C. Supernatant was aspirated and the cells were washed in 150 µl NIB and spun down again at 500G for 3 min at 4°C. This washing step was repeated once more, followed by resuspension of the cells in 2 ml Ligation Mix (400 µl 10x T4-Buffer, 40 µl 10% Tween-20, 1460 µl Annealing Buffer, and 100 µl T4 DNA Ligase, added last).

Single-cell barcoding

Using a P20 pipette, 10 µl of cell suspension in the ligation mix was added to each well of the two annealed Round 1 BC plates, taking care as not to touch the liquid at the bottom of each well. The plates were then sealed, shaken gently by hand and quickly spun down (~8 s) followed by an incubation on a shaker at 25°C for 30 min at 350 rpm. After 30 min, the cells from each well were pooled into a coated PCR strip using a P200 multichannel pipette set to 30 µl. In order to pool, each row was pipetted up and down three times before adding the liquid to a PCR strip. After eight columns were pooled into the strip, the suspension was transferred into a coated 5 ml tube on ice. This process was repeated until both plates were pooled, taking care to aspirate most liquid from the plates. The cell suspension was then spun down for 3 min at 500G at 4°C. Supernatant was aspirated and the cells were resuspended thoroughly in 2 ml new Ligation Mix. Now, 10 µl of cell suspension was added into each well of the annealed Round2 barcoding plates using a P20 pipette, taking care as to not touch the liquid within each well. The plates were sealed, shaken gently by hand and spun down quickly followed by incubating them on a shaker at 25°C for 45 min at 350 rpm. The cells were then pooled again using the above-described procedure into a new coated 15 ml tube. The cells were spun down at 500G for 3 min at 4°C. Supernatant was aspirated, the cells were washed with 150 µl NIB, and spun down again. Finally, the cells were resuspended in ~60 µl NIB (depending on total amount of cells) and counted. For counting, 5 µl of cells were mixed with 5 µl of NIB and 1x DAPI and counted on the Evos Countess II, taking the dilution into account. Sub-libraries of 3500 cells were made, and the volume was adjusted to 25 µl by addition of NIB.

Using 3500 cells results in a doublet rate of ~6.3%. The recovery rate of cells after sequencing depends on the input material (and QC thresholds), with cell lines recovering around 80% of input cells (~2800–3000 cells) and liver nuclei around 70% (~2300–2500 cells).

Reverse Crosslinking

To each sub-library of 3500 cells, 30 µl 2x Reverse Crosslinking (RC) Buffer (0.4% SDS, 100 mM NaCl, 100 mM Tris pH8.0), as well as 5 µl Proteinase K, was added. The sub-libraries were mixed and incubated on a shaker at 62°C for 1 hr at 800 rpm. Afterwards, they were transferred to a PCR cycler into a deep well module set to 62°C (lid to 80°C) for an additional hour. Afterwards, each sub-library was incubated at 80°C for 10 min and finally 5 µl of 10% Tween-20 to quench the SDS and 35 µl of NIB was added for a total volume of 100 µl.

The lysates can be stored at this point at –20°C for at least 2 days, which greatly simplifies handling many sub-libraries at once. Longer storage has not been extensively tested.

Streptavidin pull-down

Each transcript contains a biotin molecule as the RT primers are biotinylated, which is used to separate the scATAC-seq libraries from the scRNA-seq libraries. For each sub-library, 50 µl M280 Streptavidin beads were washed three times with 100 µl B&W Buffer (5 mM Tris pH 8.0, 1 M NaCl, 0.5 mM EDTA) supplemented with 0.05% Tween-20, using a magnetic stand. Afterwards, the beads were resuspended in 100 µl 2x B&W Buffer and added to the sub-library, which were then shaken at 25°C for 1 hr at 900 rpm. Now all cDNA molecules are attached to the beads whereas transposed molecules are within the supernatant. The lysate was put on a magnetic stand to separate supernatant and beads. It likely is possible to reduce the number of M280 beads in this step, significantly reducing the overall costs. However, this has not been extensively tested.

scATAC-seq library preparation

The supernatant from each sub-library was cleaned up with a Qiagen MinElute Kit and eluted twice into 30 µl 10 mM Tris pH 8.0 total. PCR Mix containing 10 µl 5x Q5 Reaction Buffer, 1 µl 10 mM dNTPs, 2 µl 10 µM i7-TruSeq-long primer, 2 µl 10 µM Nextera N5XX indexing primer, 4.5 µl H2O, and 0.5 µl Q5 polymerase was added (all oligo sequences in Supplementary file 1). Importantly, in order to distinguish the samples, each sub-library needs to be indexed with a different N5XX indexing primer. The fragments were amplified with the following protocol: 72°C for 6 min, 98°C for 1 min, then cycles of 98°C for 10 s, 66°C for 20 s, and 72°C for 45 s followed by a final incubation at 72°C for 2 min. The number of PCR cycles strongly depends on input material (liver: 17 PCR cycles, cell lines: 15 PCR cycles). The reactions were then cleaned up with custom size selection beads with 0.55x as upper cutoff and 1.4x as lower cutoff and eluted into 25 µl 10 mM Tris pH 8.0. Libraries were quantified using the Qubit HS dsDNA Quantification Kit and run on the Agilent 2100 bioanalyzer with a High Sensitivity DNA Kit.

cDNA library preparation

The beads containing the cDNA molecules were washed three times with 200 µl B&W Buffer supplemented with 0.05% Tween-20 before being resuspended in 100 µl 10 mM Tris pH 8.0 and transferred into a new PCR strip. The strip was put on a magnet and the supernatant was aspirated. The beads were then resuspended in 50 µl Template Switch Reaction Mix: 10 µl 5 X MaximaH RT Buffer, 2 µl 100 µM TS-oligo, 5 µl 10 mM dNTPs, 3 µl Enzymatics RNaseIn, 15 µl 50% PEG6000, 14 µl H2O, and 1.25 µl MaximaH RT. The sample was mixed well and incubated at 25°C for 30 min followed by an incubation at 42°C for 90 min. The beads were then washed with 100 µl 10 mM Tris while the strip was on a magnet and resuspended in 60 µl H2O. To each well, 40 µl PCR Mix was added containing 20 µl 5x Q5 Reaction Buffer, 4 µl 10 µM i7-Tru-Seq-long primer, 4 µl 10 µM Nextera N5XX Indexing primer, 2 µl 10 mM dNTPs, 9 µl H2O, and 2 µl Q5 Polymerase. The resulting mix can be split into two 50 µl PCRs or run in one 100 µl reaction. The PCR involved initial incubation at 98°C for 1 min followed by PCR cycles of 98°C for 10 s, 66°C for 20 s, and 72°C for 3 min with a final incubation at 72°C for 5 min. Importantly, in order to distinguish the samples, each sub-library needs to be indexed with a different N5XX Indexing primer. The number of PCR cycles strongly depends on input material (liver: 15 cycles, cell lines: 13 cycles).

The PCRs were cleaned up with custom size selection beads using 0.7 X as a lower cutoff (70 µl) and eluted into 25 µl 10 mM Tris pH8.0. The cDNA libraries were quantified using the Qubit HS dsDNA Quantification Kit.

scRNA-seq library preparation

As the cDNA molecules are too long for sequencing (mean length >700 bp), they need to be shortened on one side. To achieve this, 25 ng of each cDNA library was transferred to a new strip and volume was adjusted to 20 µl using H2O. Then 5 µl 5x TAPS buffer and 0.8 µl Tn5-A-only was added and the sample was incubated at 55°C for 10 min. To stop the reaction, 25 µl 1% SDS was added followed by another incubation at 55°C for 10 min. The sample was then cleaned up with custom size selection beads using a ratio of 1.3x and eluted into 30 µl. Then, 20 µl PCR mix was added containing 10 µl 5x Q5 reaction buffer, 1 µl 10 mM dNTPs, 2 µl 10 µM i7-Tru-Seq-long primer, 2 µl 10 µM Nextera N5XX Indexing primer (note: each sample needs to receive the same index primer as was used in the cDNA library preparation), 4.5 µl H2O, and 0.5 µl Q5 polymerase. The PCR was carried out with the following protocol: 72°C for 6 min, 98°C for 1 min, followed by 5 cycles of 98°C for 10 s, 66°C for 20 s, and 72°C for 45 s with a final incubation at 72°C for 2 min. Libraries were purified using custom size selection beads with a ratio of 0.5× as an upper cutoff and 0.8× as a lower cutoff. The final scRNA-seq libraries were quantified using the Qubit HS dsDNA Quantification Kit and run on the Agilent 2100 Bioanalyzer with a High Sensitivity DNA Kit.

Sequencing

Both scATAC-seq and scRNA-seq libraries were sequenced simultaneously as they were indexed with different Index 2 indices (N5XX). All libraries were sequenced on the Nova-Seq 6000 platform (Illumina) using S4 2×150 bp v1.5 kits (Read 1: 150 cycles, Index 1: 17 cycles, Index 2: 8 cycles, Read 2: 150 cycles). Libraries were partially multiplexed with standard Illumina sequencing libraries.

Custom size selection beads

To make custom size selection beads, we washed 1 ml of SpeedBeads on a magnetic stand in 1 ml of 10 mM Tris-HCl pH 8.0 and resuspended them in 50 ml Bead Buffer (9 g PEG8000, 7.3 g NaCl, 500 µl 1 M Tris HCl pH 8.0, 100 µl 0.5 M EDTA, add water to 50 ml). The beads don’t differ in their functionality from other commercially available ready-to-use size selection beads. They can be stored at 4°C for >3 months.

Analysis

Gene annotations and genomic variants

The reference genome and the Ensembl gene annotation of the C57BL/6J genome (mm10) were downloaded from Ensembl (Version GRCm38, release 102). Gene annotations for PWD/PhJ mice were downloaded from Ensembl. A consensus gene annotation set in mm10 coordinates was constructed by filtering for genes present in both gene annotations.

easySHARE-RNA-seq pre-processing

Fastq files were demultiplexed using a custom C-script, allowing one mismatch within each barcode segment (see ‘Data availability’ for the code). The reads were trimmed using cutadapt (Martin, 2011). UMIs were then extracted from bases 1–10 in Read 2 using UMI-tools (Smith et al., 2017) and added to the read name. Only reads with TTTTT at the bases 11–15 of Read 2 were kept (>96%), allowing one mismatch. Lastly, the barcode was also moved to the read name.

Species-mixing experiments

RNA-seq reads were aligned to a composite hg38-mm10 genome using STAR (Dobin et al., 2013). The resulting bamfile was then filtered for uniquely mapping reads and reads mapping to chrM, chrY, or unmapped scaffolds or containing unplaced barcodes were removed. Finally, the reads were de-duplicated using UMI-tools (Smith et al., 2017). ATAC-seq reads were also aligned to a composite genome using bwa (Li and Durbin, 2009). Duplicates were removed with Picard tools, and reads mapping to chrM, chrY, or unmapped scaffolds were filtered out. Additionally, reads that were improperly paired or had an alignment quality <30 were also removed. The reads were then split depending on which genome they mapped to, and reads per barcode were counted. Barcodes needed to be associated with at least 700 fragments and 500 UMIs in order to be considered a cell for the analysis. A barcode was considered a doublet when either the proportion of UMIs or fragments assigned to a species was less than 75%. This cutoff was chosen to mitigate possible mapping bias within the data.

easySHARE-RNA-seq processing and read alignment

We only used Read 1 for all our RNA-seq analyses as we did not need the additional genetic information for this particular analysis. Each sample was mapped to mm10 using the two-pass mode in STAR (Dobin et al., 2013) with the parameters --outFilterMultimapNmax 20 --outFilterMismatchNmax 15. We then processed the bamfiles further by moving the UMI and barcode from the read name to a bam flag, filtering out multimapping reads, and reads without a definitive barcode. To determine if a read overlapped a transcript, we used featureCounts from the subread package (Liao et al., 2014). UMI-tools was used to collapse the UMIs of aligned reads, allowing for one mismatch and de-duplication of the reads. Finally, (single-cell) count matrices were created also using UMI-tools.

easySHARE-ATAC-seq pre-processing and read alignment

Fastq files were demultiplexed using a custom C-script, allowing one mismatch within each barcode segment. The paired reads were trimmed using cutadapt (Martin, 2011), and the resulting reads were mapped to the mm10 genome using bwa mem (Li and Durbin, 2009). Reads with alignment quality <Q30, unmapped, undetermined barcode, or mapped to chrM were discarded. Duplicates were removed using Picard tools. Open chromatin regions were called by subsampling the bamfiles from all samples to a common depth, merging them into a pooled bamfile and using the peak caller MACS2 (Zhang et al., 2008) with the parameters -nomodel -keep-dup -min-length 100. The count matrices, as well as the FRiP score, were generated using featureCounts from the Subread package (Liao et al., 2014).

Filtering, integration, and dimensional reduction of scRNA-seq data

The count matrices were loaded into Seurat (Butler et al., 2018), and cells were then filtered for >200 detected genes, >500 UMIs, and <20,000 UMIs. The sub-libraries coming from the same experiment were then merged together and normalised. Merged experiments from the same species (one from male mouse, one from female mouse) were then integrated by first using SCTransform (Hafemeister and Satija, 2019) to normalise the data, then finding common features between the two experiments using FindIntegrationAnchors() and finally integrated using IntegrateData(). Lastly, the integrated datasets from C57BL/6 and PWD/PhJ were again integrated using IntegrateData(). To visualise the data, we projected the cells into 2D space by UMAP using the first 30 principal components and identified clusters using FindClusters().

Filtering, integration, and dimensional reduction of scATAC-seq data

Fragments per cell were counted using sinto, and the resulting fragment file was loaded into Signac (Stuart et al., 2021) alongside the count matrices and the peakset. We calculated basic QC statistics using base Signac, and cells were then filtered for a FRiP score of at least 0.3, >300 fragments, <15,000 fragments, a TSS enrichment >2, and a nucleosome signal <4. Again, sub-libraries coming from the same experiment were merged. We then integrated all four experiments (C57BL/6 and PWD/PhJ, one male and one female mouse each) by finding common features across datasets using FindIntegrationAnchors() using PCs 2:30 and then integrating the data using IntegrateEmbeddings(). To visualise the data, we projected the cells into 2D space by UMAP.

WNN analysis and cell-type identification

In order to use data from both modalities simultaneously, we created a multi-modal Seurat object and used WNN (Hao et al., 2021) clustering to visualise and leverage both modalities for downstream analysis. Afterwards, we assigned cell cycle scores and excluded clusters consisting of nuclei solely in the G2M-phase (2 clusters, 121 nuclei total). Cell types were assigned via expression of previously known marker genes, which allows subsetting the data into cell types. To estimate ambient RNA contamination, we used decontX (Yang et al., 2020) supplying the raw count matrix and cell-type identities. For subsequent analysis of cell-type-specific gene expression and UMAP-plots, we used the de-contaminated counts provided by decontX.

Calculating peak–gene associations

Peak–gene associations were calculated following the framework described by Ma et al., 2020. In short, Spearman correlation was calculated for every peak–gene pair within a ±500 kb window around the TSS of the expressed gene. To obtain a background estimation, we used chromVAR (Schep et al., 2017) (getBackgroundPeaks()) to generate 100 background peaks matched in GC bias and chromatin accessibility but randomly distributed throughout the genome. We calculated the Spearman correlation between every background-gene comparison, resulting in a null distribution with known population mean and standard deviation. We then calculated the z-score for the peak–gene pair in question ((correlation – population mean)/standard deviation) and used a one-sided z-test to determine the p-value. This functionality is also implemented in Signac under the function LinkPeaks(). Increasing the number of background peaks to 200, 350, or 500 for each peak–gene pair does not impact the results (data not shown).

Analysis of LSEC zonation markers

To analyse gene expression and chromatin accessibility along LSEC zonation, we ordered LSECs along pseudotime in Monocle3 (Trapnell et al., 2014). The seurat object for LSECs was converted into a Monocle cell dataset (as.cell_data_set()), clustered as a UMAP (cluster_cells()) followed by graph learning (learn_graph()). Cells with the highest Wnt2 expression were set as root for pseudotime ordering (order_cells()). Gene expression and chromatin accessibility for marker genes was smoothed over pseudotime with local polynomial regression fitting (LOESS). To identify novel marker genes, we excluded genes with low expression, divided the pseudotime into 10 bins and calculated the moving average (for three bins) across them. We then required the moving average to continuously decrease (for pericentral marker genes) or increase (for periportal marker genes), allowing two exceptions. Lastly, we divided the means for each gene by their maximum to normalise the values. Identification of CREs displaying zonation effects had equal requirements.

Gene ontology analysis

Gene ontology analysis was done using the R package clusterProfiler (Yu et al., 2012) with standard parameters.

Comparison to external datasets

In order to compare the performance of easySHARE-seq to other multiomic and single-channel technologies, we downloaded the raw data from Martin et al., 2023 (sciRNA-seq3; GEO: GSE186824, first four samples, 9.974 murine E16.5 embryonic cells), Chen et al., 2019 (SNARE-seq, GEO: GSE126074, first 4 replicates, 2.621 adult mouse brain cortex cells), Cao et al., 2018 (sciCAR, GEO:GSE117089, 5968 murine kidney cells), Cusanovich et al., 2018 (sciATAC-seq, GEO: GSE111586, 7023 murine liver nuclei), Bravo González-Blas et al., 2024 (10x Multiome ATAC+RNA; GEO: GSE218468, data generated with 10x protocol; 7.518 murine liver nuclei), and Nikopoulou et al., 2023 (10x scATAC-seq, E-MTAB: E-MTAB-12706, 7810 murine liver nuclei). We downsampled all datasets to a common sequencing depth (22,000 reads/cell for the RNA-seq, 34,000 reads/cell for the ATAC-seq) and processed them equally. To determine fragments in peaks per cell, we used the peaksets provided by the authors. Su et al., 2021a (10x 3’ expression, 82,168 murine liver cells) did not provide all raw data needed for downsampling and we therefore used the authors provided count matrix. They report a sequencing depth of ~180,000 reads/cell. Gonzales-Blas (10x Multiome, 7518 murine liver nuclei) did not report all files necessary for downsampling the ATAC-seq dataset. We therefore used the count matrix as provided. Unfortunately, Ma et al., 2020 (SHARE-seq, 42,948 murine skin cells) did neither report full raw data, code, nor sequencing depth. We therefore used the authors’ count matrix as provided.

Appendix 1

Flexibility and applicability of the easySHARE-seq framework

easySHARE-seq uses a flexible barcoding framework that can be tailored to various experimental designs. As mentioned in the main text, it allows for sequencing of fragment lengths of >200 bp, which can be critical in, e.g., studies investigating patterns of allele-specific expression or profiling of individual cancer cells and their mutations. However, in study designs not dependent on SNP coverage, sequencing costs can be cut with no downside by only sequencing 100 bp per fragment.

The entire barcoding can also be easily adapted into other protocols, such as scTCR-seq (CITR-seq; unpublished), allowing for paired investigation of T-cell receptor chains in millions of cells.

It is also straightforward to adapt easySHARE-seq to an scRNA-seq only protocol with equal or even increased throughput, as well as sample indexing, allowing to run a single experiment for, e.g., multiple replicates. To achieve this, the tagmentation step can simply be skipped and the RT primer can be switched out for /5Phos/GGGCTCGGAGATGTGTATAAGAGACAGNNNNNNNNNN-[8bp-Sample-Index]-/biodT/TTTTTTTTTTTTTTTTTTTTTTTTTVN. Before the first 10 bp UMI in R2, an 8 bp barcode was introduced which can represent an individual sample, timepoint, replicate or simply cell pool. A separate reverse transcription reaction(s) for each sample with a differently barcoded RT primer for each of them needs to be performed. Afterwards, all cells can be pooled and the protocol can be performed as described. Additionally, the amount of cells per sub-library and thus throughput can be increased linearly depending on how many different sample barcodes (RT primers) are used. For example, using two RT primers with different sample barcodes, the amount of cells per sub-library can be doubled. PCR cycles need to be adjusted for the increased cell numbers. Altogether, this allows to upscale easySHARE-seq several fold relatively straightforward with minor additional reagent costs or processing time. Throughput of easySHARE-seq is only limited by the availability of Nextera N5XX Indexing Primers, theoretically enabling the simultaneous profiling of up to a million cells.

Switching to an scATAC-seq only protocol is done by simply excluding the RT step, as well as the Streptavidin pull-down. This will significantly cut experimental cost as no RT, RNase inhibitors, or Streptavidin beads are required, bringing the cost per cell down to ~2.5 cents/cell (in a 100,000 cell experiment). Furthermore, fixation strengths can be decreased, which in our hands led to improved data quality.

However, upscaling the scATAC-seq part is less straightforward since in contrast to the RNA-seq, any additional barcodes cannot be within the insert and must be in the index reads. Currently, the barcoding scheme is designed to span the least number of bases possible which involves ‘breaking up’ the Read2 sequencing primer site between the barcoding oligos and the Tn5/RT primer. Introducing barcoded Tn5’s with additional barcodes would thus entail either (1) forgoing sample indexing on Index 2 and thus be not suitable for upscaling this part or (2) significantly increase the number of bases sequenced in Index 1 and thus decrease one of the key strengths of our protocol. Additionally, this is complicated by the need to change the sequence of the barcoding oligos, as well as the need for unassembled Tn5.

Lastly, to cut further costs on easySHARE-seq, it is possible to perform only a single ligation step. Leaving out the first ligation (in the BC plates 1) still produces easySHARE-seq libraries as the initial overhang is 8 bp long and therefore can theoretically form a stable hybridisation at room temperature.

Critical optimisation steps to use easySHARE-seq efficiently

The general molecular steps of easySHARE-seq are quite robust. However, in order to use easySHARE-seq efficiently, some prior optimisations should be performed.

As with most scRNA-seq experiments, sample preparation and fixation have the highest impact on success and quality of the experiment. As sample preparation can be quite different between tissues, general good practice is including a sufficient amount of RNase inhibitor, especially when input material is concentrated in small volumes.

The strength of fixation has a direct impact on data quality of both the scATAC-seq and scRNA-seq. In general, higher fixation leads to an increase of data quality in the scRNA-seq but makes the tagmentation in the scATAC-seq less efficient. Therefore, fixation parameters can to some extent be adjusted based on the requirements and importance of the respective output modality. Fixation strength should also be optimised in a tissue-specific manner. For example, fixing cell lines in 0.15% PFA was generally sufficient for data quality and maintaining cell integrity throughout the protocol. Liver nuclei needed a higher fixation of 0.35% and bone marrow cells (not shown) needed to be fixed in 1% PFA as they are both fragile and contain low amounts of mRNA molecules. Another critical factor is the fixation volume, e.g., fixation with 1% PFA in 1 ml leads to a different outcome than fixation in 4 ml. Generally, it is advisable to fix input material in higher volumes and with a low concentration of cells (~1 M/ml) as this leads to more consistent results and less clumping. For initial experiments to devise fixation strength, we advise to simply skip the barcoding step. This can be done by using a standard Tn5 for ATAC-seq and the following RT primer: GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG/biodT/TTTTTTTTTTTTTTTTTTTTTTTTVN. Using these modifications ensures that both the scATAC-seq and scRNA-seq can be amplified with standard Nextera N7XX and N5XX primers, allowing for cost-efficient testing of easySHARE-seq parameters. Additionally, cell integrity should be periodically checked to detect cell clumps and assess cell integrity.

Another critical aspect is minimising freeze-thaw cycles for barcoding oligos, especially for oligos containing phosphorylation modifications. Repeated freezing and thawing leads to a strong decline in protocol efficiency. Please feel free to contact the First Author for further questions.

Example workflow of easySHARE-seq library generation of 200,000 cells

To perform an experiment with a yield of ~200,000 cells, one needs to perform ~48 tagmentation reactions with 10,000 cells per reaction. After tagmentation, those get distributed into 16 RT reactions. Barcoding is then performed as described in one reaction. Afterwards, 96 sub-libraries of ~3500 cells are aliquoted and can be further processed. After Reverse Crosslinking, the samples can be stored at –20°C until the next day.

To simplify the cleanup of the lysate for the scATAC-seq library preparation, they can be cleaned up with size selection beads by adding 150 µl per well.

Discussion of ambient RNA contamination results

As seen in Figure 1—figure supplement 1, decontX identifies mean contaminated counts of 9.6% and median contaminated counts of 1.4%. The authors of decontX report mean contamination values of 1–4% in commercial droplet-based protocols and 11–14% in plate-based protocols (Yang et al., 2020), suggesting that easySHARE-seq performs better than other plate-based assays but does not reach the results of protocols that physically separate the cells. These values also suggest that few cells that are heavily contaminated strongly increase the overall estimation of contaminated counts. This could be explained by doublets and/or wrongly assigned cell types since cell-type information is passed onto decontX and builds the foundation for identifying contamination.

Overall, we found that our analyses are robust to decontamination. We thus see this filtering step as not strictly necessary but do recommend it. An alternative approach to using decontX employed in other studies is filtering out UMIs that are only associated in one sequencing read (Ma et al., 2020).

Potential alterations to the protocol in order to improve ATAC-seq data quality

Compared to the original protocol, ATAC-seq quality is decreased, potentially resulting in less resolution and variance in downstream analyses. We see multiple potential aspects of the protocol that when changed might improve ATAC-seq data quality.

The greatest potential for improvement likely lies in the fixation step as lower fixation strength generally yields higher data quality. Therefore, a simple comparison of data quality in e.g., 0.1% vs 0.2% PFA-fixated cells of the same cell type would be informative. That said, exact fixation strength is highly cell- and tissue-type dependent.

A second aspect to explore is the tagmentation step. Since fixation reduces Tn5 cutting efficiency, increasing tagmentation temperature or duration could partially compensate for this. For example, increasing tagmentation reactions to 42°C or 45 min could result in improved ATAC-seq data quality, although care must be taken not to compromise cell integrity.

Data availability

The easySHARE-seq data reported in this paper can be downloaded with the accession number GSE256434. All code used in data analysis is available at https://github.com/vosoltys/easySHARE_seq, Soltys, 2024.

The following data sets were generated
    1. Soltys V
    2. Chan YF
    (2024) NCBI Gene Expression Omnibus
    ID GSE256434. Flexible and high-throughput simultaneous profiling of gene expression and chromatin accessibility in single cells.
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE256434
The following previously published data sets were used
    1. Chen S
    2. Zhang K
    (2019) NCBI Gene Expression Omnibus
    ID GSE126074. Simultaneous profiling of transcriptome and chromatin accessibility in single nucleus.
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE126074
    1. Parekh S
    2. Tessarz P
    (2023) ArrayExpress
    ID E-MTAB-12706. scATAC-seq of mouse liver tissue derived from young and old mice.
    https://www.ebi.ac.uk/biostudies/arrayexpress/studies/E-MTAB-12706

References

    1. Forrest ARR
    2. Kawaji H
    3. Rehli M
    4. Baillie JK
    5. de Hoon MJL
    6. Haberle V
    7. Lassmann T
    8. Kulakovskiy IV
    9. Lizio M
    10. Itoh M
    11. Andersson R
    12. Mungall CJ
    13. Meehan TF
    14. Schmeier S
    15. Bertin N
    16. Jørgensen M
    17. Dimont E
    18. Arner E
    19. Schmidl C
    20. Schaefer U
    21. Medvedeva YA
    22. Plessy C
    23. Vitezic M
    24. Severin J
    25. Semple CA
    26. Ishizu Y
    27. Young RS
    28. Francescatto M
    29. Alam I
    30. Albanese D
    31. Altschuler GM
    32. Arakawa T
    33. Archer JAC
    34. Arner P
    35. Babina M
    36. Rennie S
    37. Balwierz PJ
    38. Beckhouse AG
    39. Pradhan-Bhatt S
    40. Blake JA
    41. Blumenthal A
    42. Bodega B
    43. Bonetti A
    44. Briggs J
    45. Brombacher F
    46. Burroughs AM
    47. Califano A
    48. Cannistraci CV
    49. Carbajo D
    50. Chen Y
    51. Chierici M
    52. Ciani Y
    53. Clevers HC
    54. Dalla E
    55. Davis CA
    56. Detmar M
    57. Diehl AD
    58. Dohi T
    59. Drabløs F
    60. Edge ASB
    61. Edinger M
    62. Ekwall K
    63. Endoh M
    64. Enomoto H
    65. Fagiolini M
    66. Fairbairn L
    67. Fang H
    68. Farach-Carson MC
    69. Faulkner GJ
    70. Favorov AV
    71. Fisher ME
    72. Frith MC
    73. Fujita R
    74. Fukuda S
    75. Furlanello C
    76. Furino M
    77. Furusawa J
    78. Geijtenbeek TB
    79. Gibson AP
    80. Gingeras T
    81. Goldowitz D
    82. Gough J
    83. Guhl S
    84. Guler R
    85. Gustincich S
    86. Ha TJ
    87. Hamaguchi M
    88. Hara M
    89. Harbers M
    90. Harshbarger J
    91. Hasegawa A
    92. Hasegawa Y
    93. Hashimoto T
    94. Herlyn M
    95. Hitchens KJ
    96. Ho Sui SJ
    97. Hofmann OM
    98. Hoof I
    99. Hori F
    100. Huminiecki L
    101. Iida K
    102. Ikawa T
    103. Jankovic BR
    104. Jia H
    105. Joshi A
    106. Jurman G
    107. Kaczkowski B
    108. Kai C
    109. Kaida K
    110. Kaiho A
    111. Kajiyama K
    112. Kanamori-Katayama M
    113. Kasianov AS
    114. Kasukawa T
    115. Katayama S
    116. Kato S
    117. Kawaguchi S
    118. Kawamoto H
    119. Kawamura YI
    120. Kawashima T
    121. Kempfle JS
    122. Kenna TJ
    123. Kere J
    124. Khachigian LM
    125. Kitamura T
    126. Klinken SP
    127. Knox AJ
    128. Kojima M
    129. Kojima S
    130. Kondo N
    131. Koseki H
    132. Koyasu S
    133. Krampitz S
    134. Kubosaki A
    135. Kwon AT
    136. Laros JFJ
    137. Lee W
    138. Lennartsson A
    139. Li K
    140. Lilje B
    141. Lipovich L
    142. Mackay-Sim A
    143. Manabe R
    144. Mar JC
    145. Marchand B
    146. Mathelier A
    147. Mejhert N
    148. Meynert A
    149. Mizuno Y
    150. de Lima Morais DA
    151. Morikawa H
    152. Morimoto M
    153. Moro K
    154. Motakis E
    155. Motohashi H
    156. Mummery CL
    157. Murata M
    158. Nagao-Sato S
    159. Nakachi Y
    160. Nakahara F
    161. Nakamura T
    162. Nakamura Y
    163. Nakazato K
    164. van Nimwegen E
    165. Ninomiya N
    166. Nishiyori H
    167. Noma S
    168. Noma S
    169. Noazaki T
    170. Ogishima S
    171. Ohkura N
    172. Ohimiya H
    173. Ohno H
    174. Ohshima M
    175. Okada-Hatakeyama M
    176. Okazaki Y
    177. Orlando V
    178. Ovchinnikov DA
    179. Pain A
    180. Passier R
    181. Patrikakis M
    182. Persson H
    183. Piazza S
    184. Prendergast JGD
    185. Rackham OJL
    186. Ramilowski JA
    187. Rashid M
    188. Ravasi T
    189. Rizzu P
    190. Roncador M
    191. Roy S
    192. Rye MB
    193. Saijyo E
    194. Sajantila A
    195. Saka A
    196. Sakaguchi S
    197. Sakai M
    198. Sato H
    199. Savvi S
    200. Saxena A
    201. Schneider C
    202. Schultes EA
    203. Schulze-Tanzil GG
    204. Schwegmann A
    205. Sengstag T
    206. Sheng G
    207. Shimoji H
    208. Shimoni Y
    209. Shin JW
    210. Simon C
    211. Sugiyama D
    212. Sugiyama T
    213. Suzuki M
    214. Suzuki N
    215. Swoboda RK
    216. ’t Hoen PAC
    217. Tagami M
    218. Takahashi N
    219. Takai J
    220. Tanaka H
    221. Tatsukawa H
    222. Tatum Z
    223. Thompson M
    224. Toyodo H
    225. Toyoda T
    226. Valen E
    227. van de Wetering M
    228. van den Berg LM
    229. Verado R
    230. Vijayan D
    231. Vorontsov IE
    232. Wasserman WW
    233. Watanabe S
    234. Wells CA
    235. Winteringham LN
    236. Wolvetang E
    237. Wood EJ
    238. Yamaguchi Y
    239. Yamamoto M
    240. Yoneda M
    241. Yonekura Y
    242. Yoshida S
    243. Zabierowski SE
    244. Zhang PG
    245. Zhao X
    246. Zucchelli S
    247. Summers KM
    248. Suzuki H
    249. Daub CO
    250. Kawai J
    251. Heutink P
    252. Hide W
    253. Freeman TC
    254. Lenhard B
    255. Bajic VB
    256. Taylor MS
    257. Makeev VJ
    258. Sandelin A
    259. Hume DA
    260. Carninci P
    261. Hayashizaki Y
    262. FANTOM Consortium and the RIKEN PMI and CLST (DGT)
    (2014) A promoter-level mammalian expression atlas
    Nature 507:462–470.
    https://doi.org/10.1038/nature13182

Article and author information

Author details

  1. Volker Soltys

    Friedrich Miescher Laboratory of the Max Planck Society, Tübingen, Germany
    Present address
    Max Planck Institute for Evolutionary Anthropology, Leipzig, Germany
    Contribution
    Conceptualization, Data curation, Investigation, Visualization, Writing – original draft
    For correspondence
    volker_soltys@eva.mpg.de
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0009-0008-4568-3040

  2. Moritz A Peters

    Friedrich Miescher Laboratory of the Max Planck Society, Tübingen, Germany
    Contribution
    Validation, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0009-0002-6757-3834

  3. Dingwen Su

    Friedrich Miescher Laboratory of the Max Planck Society, Tübingen, Germany
    Contribution
    Investigation, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0009-0001-4909-1279

  4. Marek Kucka

    1. Friedrich Miescher Laboratory of the Max Planck Society, Tübingen, Germany
    2. Department of Translational Genomics, University of Cologne, Cologne, Germany
    Contribution
    Investigation
    Competing interests
    No competing interests declared

  5. Yingguang Frank Chan

    1. Friedrich Miescher Laboratory of the Max Planck Society, Tübingen, Germany
    2. University of Groningen, Groningen Institute for Evolutionary Life Sciences, Groningen, Netherlands
    Contribution
    Conceptualization, Resources, Funding acquisition, Project administration, Writing – review and editing
    For correspondence
    frank.chan@rug.nl
    Competing interests
    Reviewing editor, eLife
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6292-9681

Funding

European Research Council

https://doi.org/10.3030/639096
  • Yingguang Frank Chan

European Research Council

https://doi.org/10.3030/101069216
  • Yingguang Frank Chan

Max Planck Society (International Max Planck Research School fellowship)

  • Moritz A Peters

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. Open access funding provided by Max Planck Society.

Acknowledgements

We thank members of the Chan and Jones lab for helpful discussions and critical reading of the manuscript. We are very grateful to Arnar Breevoort and Alex Pollen for sharing tissue preparation protocols and a very helpful research visit. We thank Sinja Mattes and all animal care takers at the Friedrich Miescher Laboratory for their work. We also thank the Genome Center in the Max Planck Institute for Biology Tübingen for providing support. The OP9-DL4 cells were a kind gift from Juan Carlos Zúñiga-Pflücker. MP is supported by an International Max Planck Research School fellowship. MK and YFC were supported by the European Research Council Starting Grant 639096 'HybridMiX' and Proof-of-Concept Grant 101069216 'Haplotagging'. The research was supported by the Max Planck Society.

Ethics

All animal experimental procedures were carried out under the licence number EB 01-21M at the Friedrich Miescher Laboratory of the Max Planck Society in Tübingen, Germany. The procedures were reviewed and approved by the Regierungspräsidium Tübingen, Germany.

Version history

  1. Preprint posted:
  2. Sent for peer review:
  3. Reviewed Preprint version 1:
  4. Reviewed Preprint version 2:
  5. Reviewed Preprint version 3:
  6. Version of Record published:

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.110034. This DOI represents all versions, and will always resolve to the latest one.

Copyright

© 2026, Soltys et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 843
    views
  • 71
    downloads
  • 0
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Volker Soltys
  2. Moritz A Peters
  3. Dingwen Su
  4. Marek Kucka
  5. Yingguang Frank Chan
(2026)
Flexible and high-throughput simultaneous profiling of gene expression and chromatin accessibility in single cells
eLife 15:RP110034.
https://doi.org/10.7554/eLife.110034.4

Share this article

https://doi.org/10.7554/eLife.110034

Categories and tags

Research organisms