Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

SINC-seq: correlation of transient gene expressions between nucleus and cytoplasm reflects single-cell physiology

SINC-seq: correlation of transient gene expressions between nucleus and cytoplasm reflects... We report a microfluidic system that physically separates nuclear RNA (nucRNA) and cytoplasmic RNA (cytRNA) from a single cell and enables single-cell integrated nucRNA and cytRNA-sequencing (SINC-seq). SINC-seq constructs two individual RNA-seq libraries, nucRNA and cytRNA, per cell, quantifies gene expression in the subcellular compartments, and combines them to create novel single-cell RNA-seq data. Leveraging SINC-seq, we discover distinct natures of correlation among cytRNA and nucRNA that reflect the transient physiological state of single cells. These data provide unique insights into the regulatory network of messenger RNA from the nucleus toward the cytoplasm at the single-cell level. Keywords: Single cell, RNA-seq, Microfluidics, RNA transport, Splicing, Isotachophoresis, Nucleus, Cytoplasm Background is representative of whole cells, but, to date, there has been Single-cell sequencing is a powerful tool for exploring no direct evidence of the nuclear-to-cytoplasmic correl- epigenetic, genomic, and transcriptional heterogeneities ation. Investigating the cytoplasmic RNA (cytRNA) and at an unprecedented resolution [1–6]. RNA-seq with nuclear RNA (nucRNA) expressions with the same single single nuclei (nucRNA-seq) [7, 8] is an emerging option cell are essential to assess the validity of using for profiling gene expressions of cells in tissues that nucRNA-seq, especially for the analysis of transient cannot be readily dissociated, such as the adult brain biological processes. and frozen samples. nucRNA-seq is further capable of Recent technical advances have enabled combined coupling with sorting by fluorescence-activated cell sequencing at multi-omic levels within the same single sorters [4, 7–9], Fluidigm C1 [5], and Drop-seq [10], and cells [12–14] and helped us to understand the links has demonstrated feasibilities of identifying cell types underlying the regulatory cascade. Several microfluidic and cell cycles [11] with nucRNA-seq data. Grindberg et [15–18] and non-microfluidic protocols [19, 20] offer al. [7] reported that the gene expression with single nu- parallel transcriptional and genomic analyses on the clei is similar to that with entire single cells, with only same single cell by fractionating the cytRNAs and nuclei 3.5% of the genes exhibiting differential expression. of single cells. However, we know of no work that has However, this was performed by comparing nucRNA-seq reported an integrated nucRNA-seq and cytRNA-seq vs. single-cell RNA-seq (scRNA-seq) from different single with the same cell to study RNA transport and gene cells. Such studies hypothesize that the nucRNA expression regulation and function through splicing of precursor messenger RNA (pre-mRNA) [21, 22]. Here we demonstrate a novel single-cell sequencing * Correspondence: [email protected] Mahmoud N. Abdelmoez, Kei Iida and Yusuke Oguchi contributed equally method, SINC-seq, which combines a microfluidic to this work. protocol that physically fractionates nuclear and cyto- Department of Micro Engineering, Graduate School of Engineering, Kyoto plasmic RNAs and a subcellular RNA-seq pipeline to University, Kyoto, Japan Microfluidics RIKEN Hakubi Research Team, RIKEN Cluster for Pioneering dissect RNA expressions in the individual subcellular Research, Saitama, Japan compartment. We use SINC-seq to explore both Full list of author information is available at the end of the article © The Author(s). 2018 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Abdelmoez et al. Genome Biology (2018) 19:66 Page 2 of 11 correlated and uncorrelated gene expression between per day. We hope to automate (and parallelize) the the compartments of single K562 human leukemic cells. SINC-seq protocol in thefutureand therebyincrease We further explore correlation dynamics that reflect the the rate of single-cell processing. transient response of differentiating K562 cells vs. eryth- We note that subcellular fractionation of proteins roid cells under a perturbation of sodium butyrate from single cells by electroporation was first reported by (NaB), a histone deacetylase inhibitor. We show that the Lu and co-workers [23, 24]. Our method leverages a epigenetic modification via NaB changes the degree of similar subcellular fractionation via electric field and also correlation between cytRNA and nucRNA expression uniquely enables RNA sequencing by delivering the sub- substantially. These data reveal how eukaryotes manage cellular components to two independent downstream subcellular RNA expressions via intercompartment extraction ports, including the cytRNA fraction trans- regulation. ported via ITP [16, 17]. We hope to further extend our protocol and perhaps enable protein analyses in the Results future (see Qu et al. [25] for an example of fractionation A microfluidic platform for single-cell integrated nuclear of nucleic acids vs. proteins using ITP). and cytoplasmic RNA-sequencing: SINC-seq To dissect transcriptional correlation in the subcellular Library preparation and quality control with SINC-seq compartments, we designed SINC-seq to combine elec- To critically evaluate SINC-seq, we performed experi- trophoretic fractionation of cytRNA from the nucleus ments with 93 single cells of K562 human myeloid [16–18] with off-chip RNA sequencing (Fig. 1a, b). leukemia cells and generated 186 corresponding RNA-seq SINC-seq constructs individual RNA-seq libraries with libraries using an off-chip Smart-seq2 protocol [26]. cytRNA and nucRNA and integrates the sequencing data Ziegenhain et al. [27] recently reported a comprehensive in a new format which we term an “in silico single cell.” comparison of scRNA-seq protocols including Drop-seq, SINC-seq starts with a microfluidic protocol that lever- Smart-seq with C1 (Fluidigm), and Smart-seq2. Among ages a hydrodynamic trap that captures a single cell and these methods, their work showed that Smart-seq2 is the then concentrates an electric field to selectively lyse the most sensitive with the highest number of detected genes cytoplasmic membrane while leaving the nuclear per cell. Further, Habib et al. [10, 28] recently reported a membrane relatively intact. The trap also retains the cell DroNc-seq platform approach which performs nucleus during an electric field-based extraction of single-nucleus RNA-seq. The work demonstrated that cytRNA that is initiated within 1 s of the electric field DroNc-seq detected an average of 3295 and 5134 genes, activation (see Fig. 1b, c, Additional file 1: Figure S1, respectively, for nuclei and cells of 3T3 cells. Here we have Additional file 2: Movie S1, and Methods). Hence, our leveraged the sensitivity of the Smart-seq2 protocol and a microfluidic protocol enables subcellular fractionation full-length coverage to explore the retention of introns. by coupling electric lysis [23, 24] and rapid transition to Both cytRNA-seq and nucRNA-seq of SINC-seq yielded isotachophoresis (ITP)-based nucleic acid extraction 4.64 million reads per sample (Additional file 1:Figure from single cells [16, 17]. We confirmed that the nuclei S2b, c). The average transcriptomic alignments were 94 ± retained their integrity by, for example, staining genomic 1% (mean ± standard deviation (SD)) and 93 ± 1%, respect- DNA with Hoechst (Fig. 1c). The microfluidic system ively, with cytRNA-seq and nucRNA-seq (Additional file 1: completes the entire process with voltage control via Figure S2d). Of the 93 single cells analyzed, all showed three end-channel electrodes and outputs the cytRNA successful extraction as determined by monitoring the and nucleus to different wells in less than 5 min. We ionic current of the ITP process during extraction (Add- note that the hydrodynamic trap integrated in this work itional file 1: Figure S1c). Of these 93 single cells, 84 couples hydrodynamic flow and electric field concen- passed quality control (QC) for both cytRNA-seq and tration and enables approximately 20-fold reduction nucRNA-seq. Nine of the 93 cells failed the QC for either in the applied voltage compared to our previous cytRNA-seq or nucRNA-seq. Further, in seven of the sam- protocol [16–18]. Further, the microfluidic design ples that failed QC, we observed low yield in the amplifi- allows a buffer exchange for the single cell using cation of either cytRNA or nucRNA. In two of the pressure-driven flow and introduction of the ITP buf- samples, we observed incomplete fractionation. Thus, fer chemistry. This approach helps reduce the back- after the QC, we achieved 168 data sets consisting of 84 ground noise associated with cell-free RNA in the pairs of cytRNA-seq and nucRNA-seq (see Additional file 1: original cell solution (Additional file 1:FigureS1e). Supplementary Information section titled “Fractionation These key improvements allowed us to study RNA stringency”,Additional file 1: Figure S2, Additional file 3: expressions in subcellular compartments of single TableS1, and Additionalfiles 4 and 5). cells systematically. We serially processed single cells We note that our protocol yielded smaller amounts of and prepared about eight pairs of RNA-seq libraries complementary DNA (cDNA) for extracted nucRNA Abdelmoez et al. Genome Biology (2018) 19:66 Page 3 of 11 a bc d e Fig. 1 Single-cell integrated nuclear and cytoplasmic RNA-seq (SINC-seq). a SINC-seq and conventional scRNA-seq. b Workflow of SINC-seq. Single-cell isolation at a hydrodynamic trap via pressure-driven flow (t = 0 s); lysis of cell membrane and cytRNA extraction with isotachophoresis (ITP)-aided nucleic acid extraction (t > 0 s); ITP acceleration by changing voltages (t = 40 s); voltage deactivation and sample collection from the wells of the microchannel (t >200 s). c Fluorescence microscopy images of the trapped single cell and nucleus after cytRNA extraction (stained with Hoechst) and extracted cytRNA stained with SYBR Green II. Scale bars are 20 μm. d Venn diagram of mean numbers of detected genes in cytRNA-seq and nucRNA- seq. e Percent proportion of abundance of transcripts in the cytoplasm. f Differential expression analysis between cytRNA and nucRNA. Genes enriched in cytRNA are on the right-hand side. Blue, genes with p values less than 0.001 and absolute log2 fold changes greater than unity. g Correlation coefficients of gene expression pattern computed with respect to the conventional scRNA-seq; our novel in silico single-cell normalization showed the best correlation with the scRNA-seq. We also include correlation of nucRNA vs. its in silico single cell than for cytRNA. The yield of cDNA with nucRNA was cDNA from the nucRNA fractions is due to the smaller on par with that of single nuclei prepared with an amount of RNA in a nucleus compared to the cytRNA off-the-shelf kit (PARIS Kit, Thermo Fisher Scientific) in amount for the same cell. The total amount of cDNA which the cell membrane was lysed with a chemical per single cell was 26 ± 16% less than that obtained agent. We thus hypothesize that the smaller amount of with a conventional single-cell protocol on average Abdelmoez et al. Genome Biology (2018) 19:66 Page 4 of 11 (Additional file 1: Figure S2a). We attribute this as with 0.01 < TPM < 1 in nucRNA-seq, bulk nucRNA-seq mainly due to the loss at collecting cytRNA from the detected 4034 genes. This supports the conclusion that outlet well after ITP using a standard micropipette [17]. the detection of genes was not an artifact. We further evaluated this issue by analyzing with the coverages of SINC-seq dissects the difference in subcellular gene the three lowest abundant genes (TPM~ 0.01) expression (Additional file 1:FigureS4b–d). Compared to conventional To benchmark the technical aspects of SINC-seq, we scRNA-seq, the in silico single-cell data showed a correl- assessed the sensitivity and repeatability of gene expres- ation with a mean coefficient of correlation r = 0.866 com- sion analyses with an in silico single-cell analysis. In this puted with a log-transformed expression (log10(TPM + 1)) assessment, we used 56 pairs of nucRNA-seq and (Fig. 1g). The nucRNA-seq also showed a correlation with cytRNA-seq data taken with unperturbed K562 cells scRNA-seq (different cells) with a mean coefficient of correl- which were cultured under standard conditions (without ation r = 0.711, consistent with the results of Grindberg et NaB treatment). (See a comprehensive benchmark of al. [7]. Importantly, nucRNA-seq showed a higher correl- SINC-seq in Additional file 1: Figures S3–S6 and the ation with its in silico single cell (same cells) than those with Supplementary Information section.) SINC-seq consist- single cells (different cells), supporting the conclusion that ently detected 6210 ± 1400 (mean ± SD) and 5690 ± 1500 the nucRNA-seq expression also reflects the transient gene genes per cytRNA and nucRNA, respectively, and 8200 expression of its respective single cell. This finding indicates ± 1100 genes per cell with transcripts per million (TPM) that SINC-seq for the first time reveals the correlation of greater than 1 (see Fig. 1d and comprehensive data in transient gene expression between nucRNA and cytRNA. Additional file 1: Figure S3). SINC-seq also revealed that Combined, an average gene expression profile of 12 in silico ~ 16% of transcripts were in the nucleus and ~ 84% in single cells (see Additional file 1: Figure S7 for the definition the cytoplasm (see Fig. 1e and also section titled “Com- of 12 in silico single cells) showed an excellent matching puting cyt-normalized vs. nuc-normalized data: in silico with that of 12 scRNA-seq (Pearson correlation coefficient single-cell normalization” in Methods) and also showed of r = 0.972, see Additional file 1:FigureS3q).Thetotal enriched expression of 226 and 3035 genes, respectively number of detected genes (TPM > 1) in the 12 in silico (Fig. 1f). Grindberg et al. [7] performed enrichment ana- single cells at an average sequencing depth of 2.2 M reads lyses with individually normalized reads per kilobase per was 11,131, of which 9757 genes were also detected from 12 million mapped reads (RPKM) values from nuclei and scRNA-seq averages (Additional file 1: Figure S3t). We again single cells. Compared to Grindberg’s work, our in silico stress that the in silico normalization and resulting single-cell approach uniquely enables quantitative com- scRNA-seq is a novel method which uniquely leverages a parison of gene expression and enrichment analyses physical separation and recovery as well as minimal between cytoplasm and nucleus. On average, SINC-seq cross-contamination enabled by our electrophoretic displayed about a 5.3% smaller fraction of detected fractionation. genes, 8070 ± 1100 genes at the sum total reads of 2.2 million (M) reads, than the conventional scRNA-seq, which detected 8520 ± 1100 genes (n = 12, 2.2 M reads). Cell-cycle-related genes show correlated expression in Again, we attribute this mainly to the losses associated cytoplasm and nucleus with the collection of the cytRNA from the outlet well To view the landscape of the correlation between nucRNA of the chip using a standard micropipette [17]. Notably, and cytRNA, we computed the cross-correlation of indi- our in silico single-cell data showed a wider dynamic vidual genes as a measure of covariation in the two subcel- range in the detection of genes as compared to lular compartments with 56 pairs of data taken with scRNA-seq (Additional file 1: Figure S4a). We attribute unperturbed K562 cells, and we ranked the genes by order this improvement in the dynamic range of in silico single of the value of the coefficient of correlation (Fig. 2a). We cells to the sequencing depth in nucRNA-seq. SINC-seq identified 1720 positively correlated genes and 29 nega- analyzed RNA expressions in the nucleus and cytoplasm tively correlated genes with p < 0.05. Gene ontology ana- individually with a similar number of sequencing reads. lysis revealed that the correlated genes had cell cycle as an This means that nucRNA-seq has more depth per unit enriched function (Fig. 2b) and the negative correlation RNA template than cytRNA because of the smaller had RNA splicing (Fig. 2c). We note that identification of amount of nuclear RNA. This feature facilitates the uncorrelated genes with this approach was technically detection of low abundant nuclear RNAs, hence result- challenging, as the method resulted in many uncorrelated ing in a larger dynamic range for SINC-seq in silico cell genes with negligible significance (making it impossible to analysis. To assess the validity of the gene detection with neglect the null hypothesis). In this section, we thus TPM < 1 in nucRNA-seq, we evaluated the detection of focused analyses on positively correlated and negatively genes compared with bulk nucRNA-seq. Of 4235 genes correlated genes identifiable with statistical significance. Abdelmoez et al. Genome Biology (2018) 19:66 Page 5 of 11 ab c ef g Fig. 2 Landscape of cross-correlation between cytRNA and nucRNA unveiled transcriptional oscillation of cell-cycle genes in nucRNA highly correlated with expression in cytRNA. a Quantile plot of genes sorted in order of the coefficients of cross-correlation. b Gene ontology analysis with positively correlated genes (p < 0.05) in the quantile plot and c negatively correlated genes. d–f Cell-cycle genes in in silico single-cell data, cytRNA, and nucRNA show correlation with in-phase genes (G1 vs. G1 or G2 vs. G2) and negative correlation with out-of-phase genes (G1 vs. G2). List of genes of G1 and G2 phases are provided in Additional file 1: Figure S8. g Transcriptional oscillation of cell-cycle genes in nucRNA cross- correlated with the gene expression in cytRNA To dissect the correlated gene expression, we focused To explore how the transcriptional oscillations in the on cell cycle based on transcriptional oscillations [2]and nucRNA modulated the gene expression in the cytRNA, we phase-score analysis [3](Additional file 1: Supplementary extended the analyses to compute the cross-correlation Information). Apart from the analysis of the correlation between cytRNA and nucRNA with cell-cycle genes. The landscape, we obtained a list of cell-cycle genes [29]and cell-cycle genes showed synchronized oscillation in each of extended an approach proposed by Klein et al. [2]to the two subcellular compartments (Fig. 2g), consistent with observe the behavior of individual cell-cycle genes. The in the observation that the subpopulations segregated into the silico single-cell data showed the progression of the cell G1 and G2 groups showed corresponding up-regulation cycle with a correlated variation of in-phase genes and and down-regulation of G1 and G2 genes (Additional file 1: negative correlation of out-of-phase genes (G1 vs. G2) Figure S8e, f, and Supplementary Information). Together, (Fig. 2d), consistent with scRNA-seq of K562 [2], and also these results suggest that the cytRNA and nucRNA have with the progression of the phase score (Additional file 1: similar expression patterns of cell-cycle genes and that both Figure S8a). Similarly, both cytRNA-seq and nucRNA-seq of them solely have a potency to detect the cell cycle. data revealed the cell cycle (Fig. 2e, f and Additional file 1: Figure S8b, c). Notably, we found that the negative correl- Nuclear-retained introns attenuate the transcriptional ation among out-of-phase genes was slightly higher in the oscillation − 26 nucRNA (U test, p =8×10 , Additional file 1:Figure We next studied the retained intron (RI)-mediated S8d). On the other hand, the correlation among in-phase regulation of mRNA transport [30–32] leveraging the genes was higher in the cytRNA, which may indicate intron-rich reads with nucRNA-seq of SINC-seq (see further modulation in the cytoplasm. Additional file 1: Figure S9 for comprehensive statistics Abdelmoez et al. Genome Biology (2018) 19:66 Page 6 of 11 on intron detection). After filtering RIs (Methods), genes contained three splicing-related genes [36] and a SINC-seq detected 1760 ± 570 RIs per cell, of which small nucleolar RNA (snoRNA) host gene [37]. These data 1290 ± 540 and 780 ± 290 were detected with lead us to hypothesize that the NRIs likely attenuate the nucRNA-seq and cytRNA-seq, respectively (Fig. 3a). We transcriptional oscillation in the nucleus via fine-tuning of identified 242 nuclear-retained introns (NRIs) in 213 the RNA metabolism in order to maintain the gene genes (Methods). Gene ontology analysis [33] deter- expression and functional RNAs in the cytoplasm. mined that the 213 genes had enriched functions like metabolism of RNA and RNA splicing (Fig. 3b), consist- Sodium butyrate treatment on K562 drives diverging ent with previous studies [31, 34, 35]. We examined the gene expression in subcellular compartments relationship between the probability of NRI and gene To explore the correlation dynamics under perturbation, expression in each individual (subcellular) fraction. We we further performed SINC-seq by differentiating K562 observed a positive correlation between them in the cells along the erythroid lineage with NaB (see nucRNA (r = 0.42, p < 0.01, Fig. 3c), while there was no Methods), sampling 8–13 cells per day over 5 days. We correlation in the cytRNA (r = − 0.006, p = 0.9, Fig. 3d). used NaB to modulate the histone acetylation process In contrast, we observed different enriched functions and observed transient effects of the epigenetic modifi- with cytoplasmic-enriched RIs (CRIs) and no correlation cation on the correlation through the changes in tran- between probabilities of CRI and gene expression scriptional activation and export of mRNA to the (Additional file 1: Figure S10a–c). For a better under- cytoplasm. With 41 successful SINC-seq data sets (82 standing of the function of NRIs, we examined the RNA-seq data in total), we detected differentially expression patterns of the top seven genes that were expressed genes (DEGs); 264 up-regulated and 177 highly associated with NRI-mediated regulation (Fig. 3e–g down-regulated genes were in cytRNA, and 64 and Additional file 1: Figure S10d–h). Notably, the seven up-regulated and 2 down-regulated genes were in a bc ef g Fig. 3 NRI-mediated attenuation of transcriptional oscillation in nucRNA. a Heatmap of the probability of RI in cytRNA and nucRNA fractions. NRI was identified in the upper right region indicated with the broken white line. b Gene ontology analysis with NRI, showing enriched functions like metabolism of RNA and RNA splicing. c, d Correlation analysis between the probability of NRI and the fold change of gene expression among cells with NRI and without NRI (Spl spliced) in nucRNA (upper panel) and cytRNA (lower panel), respectively. e Expression of seven genes that were highly regulated by NRI in nucRNA (p < 0.0003, Mann-Whitney U test), comparing with NRI vs. without NRI (Spl) in an individual fraction. f Coverage of SRSF5 and g HNRNPDL genes showing higher intron reads in the nuclear fraction. Coverages of ARGLU1, GAS5 (SNHG2), FBXO9, VAMP2, and PI4KAP1 genes are provided in Additional file 1: Figure S10 Abdelmoez et al. Genome Biology (2018) 19:66 Page 7 of 11 nucRNA (Additional file 1: Figure S11). We examined cytRNA and nucRNA. To leverage the different behav- the dynamics of the cross-correlation of DEG expression iors of nucRNA and cytRNA, we introduced a between cytRNA and nucRNA along the differentiation localization-embedded principal component analysis by ordering cells with the pseudotime computed using (L-PCA) that computed principal components (PCs) Monocle (version 2.4.0) [38, 39] to in silico single-cell with the subcellular gene expressions of DEGs (Fig. 4b). data (Fig. 4a, Additional file 1: Figure S12a). To provide In the L-PCA, we treated genes in cytRNA and nucRNA the statistical significance, we divided cells into five as different ones and represented a single cell as a vector groups with the pseudotime (see the divisions in Fig. 4a) having double dimension instead of using the conven- and evaluated the cross-correlation of gene expression tional approach. As expected, L-PCA resolved the trajec- between nucRNA and cytRNA (Additional file 1: Figure tory of the differentiation slightly differently compared S12b). The cross-correlation between cytRNA and to a conventional PCA that computed PCs with in silico nucRNA exhibited a gradual decrease with increasing single-cell data (Fig. 4c). To further corroborate the pseudotime, suggesting that the subcellular gene expres- L-PCA, we performed PCA on data sets of an individual sion patterns behaved differently and diverged as the fraction, showing this specific and unique transient in differentiation proceeded. We also found that non-DEGs the nucRNA — that is the day 3 cluster was furthest showed a less significant change in the coefficient of from the day 0 cluster in nucRNA (Additional file 1: Fig- cross-correlation with the differentiation (Additional file 1: ure S12f, g). With our data set, it was difficult to con- Figure S12c). We further explored this observation with clude whether L-PCA can offer better clustering than the control experiments (Additional file 1: Figure S12c–e). conventional PCA. However, these results at least dem- We found that the cross-correlation with DEGs showed onstrated that the L-PCA with SINC-seq is practical and no apparent change with unperturbed cells along the potentially useful when analyzing a biological process in- pseudotime (Additional file 1: Figure S12d). These find- volving regulations at multiple layers of single cells. ings and data with population controls at day 0 and day 4 (Additional file 1: Figure S12e) indicate that the change in Discussion the cross-correlation was not an artifact but reflects the A fundamental question is how the transcriptional oscil- dynamics of cytRNA and nucRNA along the differenti- lation in the nucRNA, which is inherently stochastic, is ation. Notably, the cross-correlation between nucRNA on transported to and correlated with gene expression in day 4 and cytRNA on day 1 (near the top right corner in the cytRNA. SINC-seq enabled direct and quantitative Fig. 4a) was lower compared to that between nucRNA on comparison of gene expressions between a nucleus, a day 1 and cytRNA on day 4 (near the bottom left corner), cytoplasm, and a whole cell of the same single cell. This suggesting that nucRNA was the driver of the diverging comparison reveals that the cell may conceivably gene expression. fine-tune a portion of its expression upon transport to We hypothesize that the divergence of gene expres- the cytoplasm (e.g., with NRI genes), while preserving sions in the two subcellular compartments reflects the correlation of other portions of its expression upon transient responses of different regulatory pathways in transport (e.g., with cell-cycle-related genes). SINC-seq nuc Pseudotime bc 01 2 3 4 L-PCA Conventional PCA (51.0) 0.8 Day4 Day2 Day1 (32.5) 10 (66.8) Day3 0.6 Day0 Day2 Day0 0.4 Day1 −10 Day3 Day4 (67.4) −20 0.2 −20 −30 −25 0 25 −25 0 25 PC1 PC1 Fig. 4 Differentiation of K562 cells to erythroid cells shows a dynamical change of cross-correlation between cytRNA and nucRNA. a The cross- correlation of DEG expression between cytRNA and nucRNA along with the pseudotime. Numbers along axes show groups of cells used in the analysis shown in Additional file 1: Figure S12b. b L-PCA analysis of differentiating K562 cells to erythroid cells by sodium butyrate treatment. c Conventional PCA analysis with differentiating K562 cells cyt Pseudotime 43 2 1 0 PC2 PC2 Abdelmoez et al. Genome Biology (2018) 19:66 Page 8 of 11 also revealed that the cells under external perturbation Buffers dynamically alter the correlation and exhibit a unique We designed buffers for ITP-based selective extraction, trajectory of differentiation at subcellular resolution. separation (from the trapped nucleus), purification, and These findings shed new light on the characteristics of transport of cytRNA to the cytRNA output well of the post-transcriptional regulations with a single cell and chip (see more details in Shintaku et al. [16] and subcellular compartment resolution. Kuriyama et al. [17]). The leading electrolyte buffer (LE) Our study also suggests a compelling caution to an components were 50 mM Tris and 25 mM HCl contain- approach that approximates the transcriptomic profile of ing 0.4% poly(vinylpyrrolidone) (PVP) (calculated pH of the whole cell with that of a single compartment without 8.1). The trailing electrolyte buffer (TE) components validation. The SINC-seq platform will be broadly were 50 mM imidazole and 25 mM HEPES containing applicable to different types of cells as long as they are (initial calculated pH of 8.3) 0.4% PVP. We included isolated as singles. The method thus will contribute to PVP to suppress electroosmotic flow. We purchased validate existing subcellular RNA-seq methods [4, 5, 9– Tris, HEPES, imidazole, and HCl from Sigma-Aldrich, 11, 19, 20] and also define their limitations. and PVP (molecular weight 1 MDa) from Polyscience. We prepared all solutions in UltraPure DNase-/RNase-free deionized (DI) water (Life Technologies). Conclusion To dissect transcriptional correlation in subcellular com- Microfluidic system setup partments, we devised SINC-seq, which enables inte- We fabricated polydimethylsiloxane (PDMS, Sylgard grated nuclear and cytoplasmic RNA-seq of single cells 184, Dow Corning) microchannel superstructures by coupling physical fractionation of cytRNA from the (Additional file 1: Figure S1a, b) with a soft lithography nucleus of a single cell with a high-throughput and bonded them to a glass substrate. SU-8 (SU-8 2025, RNA-seq. Leveraging SINC-seq, we explored the land- MicroChem) molds were prepared on glass substrates scape of the correlation between nucRNA and cytRNA with the microchannel patterns made of chromium thin with a total of 84 K562 cells, which corresponds to 168 films, exposing the SU-8 to ultraviolet light through the RNA-seq libraries. The SINC-seq data unveiled three pattern. The nominal channel width and depth of the distinct natures of correlation among cytRNA and microchannels were 50 μm and 25 μm, respectively. We nucRNA that reflected the physiological state of single designed 3-μm-wide and 5-μm-long hydrodynamic traps. cells: highly correlated expression in cell-cycle-related We have optimized the size of the hydrodynamic trap as genes, the distorted correlation via NRIs, and the correl- 3 μm wide and 25 μm deep so that it can capture a ation dynamics along the differentiation of K562 cells to single cell and trap a nucleus during the electric field erythroid cells under sodium butyrate perturbation. extraction of cytRNA. We observed some deformation These data provide unique insights into the regulatory of cells due to the presence of the pressure-driven flow, network of mRNA from the nucleus toward the cyto- but observed no mechanical lyses prior to the applica- plasm at the single-cell level. tion of the electric field. This trapping process had a timescale of several minutes. Methods Before each experiment, we preconditioned the micro- Cells channel by filling the inlet and outlet wells with washing We purchased K562 cells (human lymphoblast, solutions and applying a vacuum to the waste well. Our chronic myelogenous leukemia) from RIKEN BioRe- washing process was as follows: 1 M NaOH for 1 min, source Center and the Japanese Collection of 1 M HCl for 1 min, and DI water for 1 min. All washing Research Bioresources (JCRB) cell bank. We cultured solutions contained 0.1% Triton X-100 to suppress the K562 cells in RPMI-1640 Medium (Life Technolo- bubble clogging in the hydrophobic microchannel. gies) with 10% fetal bovine serum and 1% Following this, we loaded 9.5 μL of LE and TE to the penicillin-streptomycin-glutamine at 37 °C in 5% CO . outlet and inlet wells, respectively, and briefly applied a We washed the cells with phosphate-buffered saline vacuum to the waste well to exchange the solution in once and suspended them in a sample buffer the microchannel with LE and TE. The hydrodynamic containing 50 mM imidazole, 25 mM 4-(2-hydro- pressure induced by buffers in the inlet and outlet wells xyethyl)-1-piperazineethanesulfonic acid (HEPES), and created a pressure-driven laminar flow from both inlet 175 mM sucrose (pH 7.6) at a concentration of ~ 0.8 and outlet wells toward the waste well and formed a cells/μL and stored them on ice until the experiments stable LE-TE interface at the junction of three micro- were performed. To differentiate the K562 cells, we channels. We then picked up a single cell of interest incubated them with 1 mM NaB (Sigma-Aldrich, from the cell suspension by aspirating 1 μL using a B5887) and harvested them after 96 h of induction. standard micropipette. We observed this process using a Abdelmoez et al. Genome Biology (2018) 19:66 Page 9 of 11 microscope and confirmed aspiration of a single cell. using a Nextera XT DNA Sample Preparation Kit (Illu- We then dispensed this 1-μL volume into the inlet well mina) and purified the cDNA library following the man- and introduced it into the microchannel via the ufacturer’s protocol, except that we eluted the cDNA pressure-driven flow. Once we visually confirmed the sample with 24 μL instead of 50 μL (see Additional file 1: captured single cell at the hydrodynamic trap (Fig. 1c), Figure S2a for yields of cDNA). We pooled 98–108 we added 9.5 μL of the TE to the waste well to reduce libraries and sequenced them on an Illumina HiSeq the pressure-driven flow. We aborted our protocol in 2500 with 100-base paired-end reads to an average depth cases where we observed two or more cells at the hydro- of 4.64 million reads (Additional file 1: Figure S2b, c). dynamic trap. We placed 300-μm-diameter platinum We mapped the trimmed sequencing reads to the tran- wire electrodes into the wells and applied − 150 V, − scripts derived from the human reference genome 170 V, and 0 V to the electrodes at the inlet, waste, and (GRCh37.75) using the STAR (version 2.5.1b) mapping outlet wells, respectively. The DC voltage created a con- program [40] with ENCODE options, and calculated centrated electric field at the hydrodynamic trap expression estimates with TPM using RNA-Seq by (Additional file 1: Figure S1d) and lysed the cytoplasmic Expectation-Maximization (RSEM v1.3.0) [41]. membrane within 1 s. Appropriate placement of the ITP We performed DE analyses individually with cytRNA buffers with the DC electric field enabled an immediate and nucRNA comparing day 1–day 4 samples to day 0 transition from the lysis to an ITP process that collects samples. We used the MATLAB functions “nbintest” and focuses cytRNA into an ITP-zone, TE-to-LE inter- and “mafdr” and identified significance with p values less face. At 40 s, we changed the voltages to − 350 V and − than 0.001 and absolute log2 fold changes greater than 510 V at the inlet and waste wells, respectively, to accel- unity. erate the migration of the ITP zone. The ITP zone trans- ported the cytRNA to the output well in about 100 s, Analysis of intron retention while the nucleus was retained at the hydrodynamic We computed intron expressions with fragments per trap. We also monitored the current during the extrac- kilobase of intron per million mapped reads (FPKM) tion with a computer running a custom MATLAB using 347,041 unique introns (longer than 50 nt) on the (Mathworks, Inc., Natick, MA, USA) script. The magni- genome annotation with 56 SINC-seq data of K562 cells tude of the current decreased as the ITP zone (contain- under standard culturing conditions. On average, ing the focused cytRNA) advanced in the channel and as SINC-seq yielded 14.8% reads mapped to introns with the lower conductivity TE replaced the higher conduct- nucRNA-seq, but only 1.1% with cytRNA-seq. Further, ivity LE (Additional file 1: Figure S1c). The current SINC-seq detected 38,800 ± 10,000 and 34,800 ± 8000 signal plateaued near t = 100 s, coincident with the time unique introns in nucRNA-seq and cytRNA-seq, re- at which the focused cytRNA eluted into the outlet well. spectively, and 56,300 ± 9500 per cell with FPKM of We deactivated the voltages at 200 s and used a standard more than 0. We identified an RI that had at least 10% pipette to transfer two aliquots from the chip: 9.5 μL expression level of the gene, 95% coverage in the in- from the outlet well containing the cytRNA and 1 μL tronic region, and non-zero expression in the adjacent containing the cell nucleus from the inlet well. Detailed exon. We discarded intron reads locating on a gene with descriptions of a similar protocol and chip, together with less than 2 TPM. On the other hand, we identified a a narrated video description, were reported by Kuriyama fully spliced intron that had less than 1% expression of et al. [18]. the gene and 50% coverage in the adjacent exon. We dis- carded intron reads that failed the preceding criteria. Library preparation and mapping analysis We validated the RI identification with splice site scores We synthesized respective cDNA libraries from the [42], which showed lower values with RIs than fully fractionated cytRNA and nucRNA separately using spliced introns, using 9mer (exonic 3mer + intronic − 16 Smart-seq2 (SMART-seq v4 Ultra Low Input RNA Kit 6mer) around the 5′ splice site (p value < 2.2 × 10 , U for Sequencing, Clontech) with 18 polymerase chain test), and 23mer (intronic 20mer + exonic 3mer) around − 12 reaction (PCR) cycles followed by purification with the 3′ splice site (p value < 2.5 × 10 , U test). In total, Agencourt AMPure XP (Beckman Coulter). We exam- we detected 17,277 RIs, of which 14,134 and 2316 RIs ined the yield and quality of cDNA, respectively, with a had a higher probability of RI in nucRNA and cytRNA, Qubit 2.0 Fluorometer (Thermo Fisher Scientific) and respectively. The RI enrichment in the nucleus was quantitative PCR (qPCR) targeting GAPDH (glyceralde- consistent with that for previous studies [30–32]. hyde-3-phosphate dehydrogenase, Hs02758991_g1, To identify NRIs, we calculated the probability of in- Thermo Fisher Scientific) and HBG (gamma-globin tron retention defined as the proportion of cells with the genes, Hs00361131_g1, Thermo Fisher Scientific). We RI and identified NRIs that had 0.25 higher probability performed the tagmentation reaction with 200 pg cDNA in the nuclear fraction than in the cytoplasmic fraction. Abdelmoez et al. Genome Biology (2018) 19:66 Page 10 of 11 We then filtered unique NRIs, discarding smaller NRIs Additional files that had overlap with a long NRI. On the other hand, we Additional file 1: Supplementary information and Figures S1–S12. identified CRIs that had 0.25 higher probability in the (DOCX 4587 kb) cytoplasmic fraction than in the nuclear fraction. Additional file 2: Movie S1. Electrical lysis and RNA extraction visualized by SYBR Green II. (MOV 1279 kb) Computing cyt-normalized vs. nuc-normalized data: in Additional file 3: Table S1. Quality control of SINC-seq samples. silico single-cell normalization (XLSX 21 kb) We computed in silico single-cell RNA-seq (scRNA-seq) Additional file 4: Movie S2. Experimental run #55 with unsuccessful fractionation. (MP4 1372 kb) data with cytRNA-seq and nucRNA-seq data, scaling the Additional file 5: Movie S3. Experimental run #69 with unsuccessful raw TPM values and combining the cytRNA-seq and fractionation. (MP4 1116 kb) nucRNA-seq data as Acknowledgements TPM ¼ TPM þ TPM ; ðS1Þ in−silico cyt nuc We acknowledge Yuka Ozaki of Kyoto University for technical support in the preparation of the sequencing libraries. TPM ¼ αTPM ; ðS2Þ cyt Funding This work was funded by the ImPACT Program of the Council for Science, TPM ¼ βTPM ; ðS3Þ nuc Technology, and Innovation (Cabinet Office, Government of Japan) and also by the Japan Society for the Promotion of Science under 26289035 and where α and β are normalization factors, which make the summation of the TPM values, TPM , of the in in silico Availability of data and materials silico single-cell data to be one million. We here write The data from this study have been deposited into the Sequencing Read Archive and are available for download under accession SRP119800 [43]. raw TPM values with an asterisk and TPM values of in silico single-cell data, cytRNA-seq, and nucRNA-seq, Authors’ contributions respectively, with their subscripts. We calculate α and β MNA performed the experiments; YO, KI, HS, and HN analyzed the data. SU, RY, and HK helped to analyze the data. HS and JGS designed and supervised using ΔCt of qPCR data taken at the QC of cDNA the project. All authors wrote the manuscript. All authors read and approved (Additional file 3: Table S1) as the final manuscript. ΔCt 2 Ethics approval α ¼ ; ðS4Þ Does not apply to this study. TPM cyt;geneA ΔCt þ 2 TPM nuc;geneA Competing interests H.S., Y.O., and S.U. submitted patent applications (Japanese Patent Application No. 2016-014708, No. 2016-177,163, and PCT/JP2017/003069) β ¼ 1−α; ðS5Þ relating to this work. Correspondence and requests for materials should be addressed to H.S. ([email protected]). * * where TPM , TPM , and ΔCt are, respect- cyt,geneA nuc, geneA ively, the raw TPM value of gene A with cytRNA-seq, Publisher’sNote Springer Nature remains neutral with regard to jurisdictional claims in the raw TPM value of gene A with nucRNA-seq, and published maps and institutional affiliations. ΔCt = Ct – Ct which is ΔCt with nuc, geneA cyt, geneA, respect to gene A. We calculated pairs of α and β with Author details Department of Micro Engineering, Graduate School of Engineering, Kyoto GAPDH and HBG genes (HBG1 + HBG2)as gene As, University, Kyoto, Japan. Microfluidics RIKEN Hakubi Research Team, RIKEN and used the mean α and β to compute the TPM. Cluster for Pioneering Research, Saitama, Japan. Medical Research Support We here note that α and β, respectively, indicate the Center, Graduate School of Medicine, Kyoto University, Kyoto, Japan. Department of Biological Sciences, Graduate School of Science, The (complementary) fractions of cytoplasmic transcripts University of Tokyo, Tokyo, Japan. Department of Hematology, Faculty of and nuclear transcripts in the in silico single cells. We Medicine, University of Tsukuba, Tsukuba, Japan. Department of Mechanical thus can quantify the fraction of the cytoplasmic tran- Engineering, Stanford University, Stanford, CA, USA. script as shown in Fig. 1e. Received: 13 October 2017 Accepted: 7 May 2018 L-PCA The PCA for the conventional scRNA-seq uses an (n × m) References 1. Wu AR, et al. Quantitative assessment of single-cell RNA-sequencing matrix of gene expressions (n = number of genes) with methods. Nat Methods. 2014;11(1):41–6. multiple samples (m = number of samples); however, our 2. Klein AM, et al. Droplet barcoding for single-cell transcriptomics applied to L-PCA uses a (2n × m)matrix of geneexpressions,having embryonic stem cells. Cell. 2015;161(5):1187–201. 3. Macosko EZ, et al. Highly parallel genome-wide expression profiling of a twofold dimension compared to the conventional individual cells using nanoliter droplets. Cell. 2015;161(5):1202–14. matrix, derived from cytRNA-seq and nucRNA-seq. The 4. Habib N, et al. Div-Seq: single-nucleus RNA-Seq reveals dynamics of rare L-PCA was performed using PCA with “prcomp” in R. adult newborn neurons. Science. 2016;353(6302):925–8. Abdelmoez et al. Genome Biology (2018) 19:66 Page 11 of 11 5. Lake BB, et al. Neuronal subtypes and diversity revealed by single-nucleus 37. Williams GT, Farzaneh F. Are snoRNAs and snoRNA host genes new players RNA sequencing of the human brain. Science. 2016;352(6293):1586–90. in cancer? Nat Rev Cancer. 2012;12(2):84–8. 6. Tanay A, Regev A. Scaling single-cell genomics from phenomenology to 38. Qiu X, et al. Single-cell mRNA quantification and differential analysis with mechanism. Nature. 2017;541(7637):331–8. Census. Nat Methods. 2017;14(3):309–15. 7. Grindberg RV, et al. RNA-sequencing from single nuclei. Proc Natl Acad Sci 39. Trapnell C, et al. The dynamics and regulators of cell fate decisions are U S A. 2013;110(49):19802–7. revealed by pseudotemporal ordering of single cells. Nat Biotechnol. 2014; 32(4):381–6. 8. Krishnaswami SR, et al. Using single nuclei for RNA-seq to capture the 40. Dobin A, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. transcriptome of postmortem neurons. Nat Protoc. 2016;11(3):499–524. 2013;29(1):15–21. 9. Lacar B, et al. Nuclear RNA-seq of single neurons reveals molecular 41. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data signatures of activation. Nat Commun. 2016;7:11022. with or without a reference genome. BMC Bioinformatics. 2011;12:323. 10. Habib N, et al. Massively parallel single-nucleus RNA-seq with DroNc-seq. 42. Yeo G, Burge CB. Maximum entropy modeling of short sequence motifs with Nat Methods. 2017;14(10):955–8. applications to RNA splicing signals. J Comput Biol. 2004;11(2-3):377–94. 11. Gao R, et al. Nanogrid single-nucleus RNA sequencing reveals phenotypic 43. Abdelmoez MN, et al., Gene expressions in nucleus and cytoplasm at the diversity in breast cancer. Nat Commun. 2017;8(1):228. single cell resolution. Sequencing Read Archive. 2017. https://www.ncbi.nlm. 12. Macaulay IC, et al. G&T-seq: parallel sequencing of single-cell genomes and nih.gov/bioproject/PRJNA413830. transcriptomes. Nat Methods. 2015;12(6):519–22. 13. Angermueller C, et al. Parallel single-cell sequencing links transcriptional and epigenetic heterogeneity. Nat Methods. 2016;13(3):229–32. 14. Cheow LF, et al. Single-cell multimodal profiling reveals cellular epigenetic heterogeneity. Nat Methods. 2016;13(10):833–6. 15. Han L, et al. Co-detection and sequencing of genes and transcripts from the same single cells facilitated by a microfluidics platform. Sci Rep. 2014;4:6485. 16. Shintaku H, et al. On-chip separation and analysis of RNA and DNA from single cells. Anal Chem. 2014;86(4):1953–7. 17. Kuriyama K, Shintaku H, Santiago JG. Isotachophoresis for fractionation and recovery of cytoplasmic RNA and nucleus from single cells. Electrophoresis. 2015;36(14):1658–62. 18. Kuriyama K, Shintaku H, Santiago JG. Protocol for microfluidic system to automate the preparation and fractionation of the nucleic acids in the cytoplasm versus nuclei of single cells. Bio-protocol. 2016;6(12):e1844. 19. Hou Y, et al. Single-cell triple omics sequencing reveals genetic, epigenetic, and transcriptomic heterogeneity in hepatocellular carcinomas. Cell Res. 2016;26(3):304–19. 20. Hu Y, et al. Simultaneous profiling of transcriptome and DNA methylome from a single cell. Genome Biol. 2016;17:88. 21. Kohler A, Hurt E. Exporting RNA from the nucleus to the cytoplasm. Nat Rev Mol Cell Biol. 2007;8(10):761–73. 22. Gerstberger S, Hafner M, Tuschl T. A census of human RNA-binding proteins. Nat Rev Genet. 2014;15(12):829–45. 23. Wang J, Lu C, et al. Kinetics of NF-κB nucleocytoplasmic transport probed by single-cell screening without imaging. Lab Chip. 2010;10(21):2911–6. 24. Zhan Y, Lu C, et al. One-step extraction of subcellular proteins from eukaryotic cells. Lab Chip. 2010;10(16):2046–8. 25. Qu, Marshall LA, Santiago JG. Simultaneous purification and fractionation of nucleic acids and proteins from complex samples using bidirectional isotachophoresis. Anal Chem. 2014;86(15):7264–8. 26. Picelli S, et al. Smart-seq2 for sensitive full-length transcriptome profiling in single cells. Nat Methods. 2013;10(11):1096–8. 27. Ziegenhain C, et al. Comparative analysis of single-cell RNA sequencing methods. Mol Cell. 2017;65(4):631–43. e4 28. Habib N, et al. DroNc-Seq: deciphering cell types in human archived brain tissues by massively-parallel single nucleus RNA-seq. bioRxiv. 2017; https:// doi.org/10.1101/115196 29. Whitfield ML, et al. Identification of genes periodically expressed in the human cell cycle and their expression in tumors. Mol Biol Cell. 2002;13(6):1977–2000. 30. Ninomiya K, Kataoka N, Hagiwara M. Stress-responsive maturation of Clk1/4 pre-mRNAs promotes phosphorylation of SR splicing factor. J Cell Biol. 2011;195(1):27–40. 31. Boutz PL, Bhutkar A, Sharp PA. Detained introns are a novel, widespread class of post-transcriptionally spliced introns. Genes Dev. 2015;29(1):63–80. 32. Naro C, et al. An orchestrated intron retention program in meiosis controls timely usage of transcripts during germ cell differentiation. Dev Cell. 2017; 41(1):82–93. e4 33. Tripathi S, et al. Meta- and orthogonal integration of influenza “OMICs” data defines a role for UBR4 in virus budding. Cell Host Microbe. 2015;18(6):723–35. 34. Iida K, Go M. Survey of conserved alternative splicing events of mRNAs encoding SR proteins in land plants. Mol Biol Evol. 2006;23(5):1085–94. 35. Hsu TY, et al. The spliceosome is a therapeutic vulnerability in MYC-driven cancer. Nature. 2015;525(7569):384–8. 36. Lareau LF, et al. Unproductive splicing of SR genes associated with highly conserved and ultraconserved DNA elements. Nature. 2007;446(7138):926–9. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Genome Biology Springer Journals

SINC-seq: correlation of transient gene expressions between nucleus and cytoplasm reflects single-cell physiology

Download PDF
11 pages

Loading next page...
 
/lp/springer_journal/sinc-seq-correlation-of-transient-gene-expressions-between-nucleus-and-gtExPkSNde

References (46)

Publisher
Springer Journals
Copyright
Copyright © 2018 by The Author(s).
Subject
Life Sciences; Animal Genetics and Genomics; Human Genetics; Plant Genetics and Genomics; Microbial Genetics and Genomics; Bioinformatics; Evolutionary Biology
eISSN
1474-760X
DOI
10.1186/s13059-018-1446-9
pmid
29871653
Publisher site
See Article on Publisher Site

Abstract

We report a microfluidic system that physically separates nuclear RNA (nucRNA) and cytoplasmic RNA (cytRNA) from a single cell and enables single-cell integrated nucRNA and cytRNA-sequencing (SINC-seq). SINC-seq constructs two individual RNA-seq libraries, nucRNA and cytRNA, per cell, quantifies gene expression in the subcellular compartments, and combines them to create novel single-cell RNA-seq data. Leveraging SINC-seq, we discover distinct natures of correlation among cytRNA and nucRNA that reflect the transient physiological state of single cells. These data provide unique insights into the regulatory network of messenger RNA from the nucleus toward the cytoplasm at the single-cell level. Keywords: Single cell, RNA-seq, Microfluidics, RNA transport, Splicing, Isotachophoresis, Nucleus, Cytoplasm Background is representative of whole cells, but, to date, there has been Single-cell sequencing is a powerful tool for exploring no direct evidence of the nuclear-to-cytoplasmic correl- epigenetic, genomic, and transcriptional heterogeneities ation. Investigating the cytoplasmic RNA (cytRNA) and at an unprecedented resolution [1–6]. RNA-seq with nuclear RNA (nucRNA) expressions with the same single single nuclei (nucRNA-seq) [7, 8] is an emerging option cell are essential to assess the validity of using for profiling gene expressions of cells in tissues that nucRNA-seq, especially for the analysis of transient cannot be readily dissociated, such as the adult brain biological processes. and frozen samples. nucRNA-seq is further capable of Recent technical advances have enabled combined coupling with sorting by fluorescence-activated cell sequencing at multi-omic levels within the same single sorters [4, 7–9], Fluidigm C1 [5], and Drop-seq [10], and cells [12–14] and helped us to understand the links has demonstrated feasibilities of identifying cell types underlying the regulatory cascade. Several microfluidic and cell cycles [11] with nucRNA-seq data. Grindberg et [15–18] and non-microfluidic protocols [19, 20] offer al. [7] reported that the gene expression with single nu- parallel transcriptional and genomic analyses on the clei is similar to that with entire single cells, with only same single cell by fractionating the cytRNAs and nuclei 3.5% of the genes exhibiting differential expression. of single cells. However, we know of no work that has However, this was performed by comparing nucRNA-seq reported an integrated nucRNA-seq and cytRNA-seq vs. single-cell RNA-seq (scRNA-seq) from different single with the same cell to study RNA transport and gene cells. Such studies hypothesize that the nucRNA expression regulation and function through splicing of precursor messenger RNA (pre-mRNA) [21, 22]. Here we demonstrate a novel single-cell sequencing * Correspondence: [email protected] Mahmoud N. Abdelmoez, Kei Iida and Yusuke Oguchi contributed equally method, SINC-seq, which combines a microfluidic to this work. protocol that physically fractionates nuclear and cyto- Department of Micro Engineering, Graduate School of Engineering, Kyoto plasmic RNAs and a subcellular RNA-seq pipeline to University, Kyoto, Japan Microfluidics RIKEN Hakubi Research Team, RIKEN Cluster for Pioneering dissect RNA expressions in the individual subcellular Research, Saitama, Japan compartment. We use SINC-seq to explore both Full list of author information is available at the end of the article © The Author(s). 2018 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Abdelmoez et al. Genome Biology (2018) 19:66 Page 2 of 11 correlated and uncorrelated gene expression between per day. We hope to automate (and parallelize) the the compartments of single K562 human leukemic cells. SINC-seq protocol in thefutureand therebyincrease We further explore correlation dynamics that reflect the the rate of single-cell processing. transient response of differentiating K562 cells vs. eryth- We note that subcellular fractionation of proteins roid cells under a perturbation of sodium butyrate from single cells by electroporation was first reported by (NaB), a histone deacetylase inhibitor. We show that the Lu and co-workers [23, 24]. Our method leverages a epigenetic modification via NaB changes the degree of similar subcellular fractionation via electric field and also correlation between cytRNA and nucRNA expression uniquely enables RNA sequencing by delivering the sub- substantially. These data reveal how eukaryotes manage cellular components to two independent downstream subcellular RNA expressions via intercompartment extraction ports, including the cytRNA fraction trans- regulation. ported via ITP [16, 17]. We hope to further extend our protocol and perhaps enable protein analyses in the Results future (see Qu et al. [25] for an example of fractionation A microfluidic platform for single-cell integrated nuclear of nucleic acids vs. proteins using ITP). and cytoplasmic RNA-sequencing: SINC-seq To dissect transcriptional correlation in the subcellular Library preparation and quality control with SINC-seq compartments, we designed SINC-seq to combine elec- To critically evaluate SINC-seq, we performed experi- trophoretic fractionation of cytRNA from the nucleus ments with 93 single cells of K562 human myeloid [16–18] with off-chip RNA sequencing (Fig. 1a, b). leukemia cells and generated 186 corresponding RNA-seq SINC-seq constructs individual RNA-seq libraries with libraries using an off-chip Smart-seq2 protocol [26]. cytRNA and nucRNA and integrates the sequencing data Ziegenhain et al. [27] recently reported a comprehensive in a new format which we term an “in silico single cell.” comparison of scRNA-seq protocols including Drop-seq, SINC-seq starts with a microfluidic protocol that lever- Smart-seq with C1 (Fluidigm), and Smart-seq2. Among ages a hydrodynamic trap that captures a single cell and these methods, their work showed that Smart-seq2 is the then concentrates an electric field to selectively lyse the most sensitive with the highest number of detected genes cytoplasmic membrane while leaving the nuclear per cell. Further, Habib et al. [10, 28] recently reported a membrane relatively intact. The trap also retains the cell DroNc-seq platform approach which performs nucleus during an electric field-based extraction of single-nucleus RNA-seq. The work demonstrated that cytRNA that is initiated within 1 s of the electric field DroNc-seq detected an average of 3295 and 5134 genes, activation (see Fig. 1b, c, Additional file 1: Figure S1, respectively, for nuclei and cells of 3T3 cells. Here we have Additional file 2: Movie S1, and Methods). Hence, our leveraged the sensitivity of the Smart-seq2 protocol and a microfluidic protocol enables subcellular fractionation full-length coverage to explore the retention of introns. by coupling electric lysis [23, 24] and rapid transition to Both cytRNA-seq and nucRNA-seq of SINC-seq yielded isotachophoresis (ITP)-based nucleic acid extraction 4.64 million reads per sample (Additional file 1:Figure from single cells [16, 17]. We confirmed that the nuclei S2b, c). The average transcriptomic alignments were 94 ± retained their integrity by, for example, staining genomic 1% (mean ± standard deviation (SD)) and 93 ± 1%, respect- DNA with Hoechst (Fig. 1c). The microfluidic system ively, with cytRNA-seq and nucRNA-seq (Additional file 1: completes the entire process with voltage control via Figure S2d). Of the 93 single cells analyzed, all showed three end-channel electrodes and outputs the cytRNA successful extraction as determined by monitoring the and nucleus to different wells in less than 5 min. We ionic current of the ITP process during extraction (Add- note that the hydrodynamic trap integrated in this work itional file 1: Figure S1c). Of these 93 single cells, 84 couples hydrodynamic flow and electric field concen- passed quality control (QC) for both cytRNA-seq and tration and enables approximately 20-fold reduction nucRNA-seq. Nine of the 93 cells failed the QC for either in the applied voltage compared to our previous cytRNA-seq or nucRNA-seq. Further, in seven of the sam- protocol [16–18]. Further, the microfluidic design ples that failed QC, we observed low yield in the amplifi- allows a buffer exchange for the single cell using cation of either cytRNA or nucRNA. In two of the pressure-driven flow and introduction of the ITP buf- samples, we observed incomplete fractionation. Thus, fer chemistry. This approach helps reduce the back- after the QC, we achieved 168 data sets consisting of 84 ground noise associated with cell-free RNA in the pairs of cytRNA-seq and nucRNA-seq (see Additional file 1: original cell solution (Additional file 1:FigureS1e). Supplementary Information section titled “Fractionation These key improvements allowed us to study RNA stringency”,Additional file 1: Figure S2, Additional file 3: expressions in subcellular compartments of single TableS1, and Additionalfiles 4 and 5). cells systematically. We serially processed single cells We note that our protocol yielded smaller amounts of and prepared about eight pairs of RNA-seq libraries complementary DNA (cDNA) for extracted nucRNA Abdelmoez et al. Genome Biology (2018) 19:66 Page 3 of 11 a bc d e Fig. 1 Single-cell integrated nuclear and cytoplasmic RNA-seq (SINC-seq). a SINC-seq and conventional scRNA-seq. b Workflow of SINC-seq. Single-cell isolation at a hydrodynamic trap via pressure-driven flow (t = 0 s); lysis of cell membrane and cytRNA extraction with isotachophoresis (ITP)-aided nucleic acid extraction (t > 0 s); ITP acceleration by changing voltages (t = 40 s); voltage deactivation and sample collection from the wells of the microchannel (t >200 s). c Fluorescence microscopy images of the trapped single cell and nucleus after cytRNA extraction (stained with Hoechst) and extracted cytRNA stained with SYBR Green II. Scale bars are 20 μm. d Venn diagram of mean numbers of detected genes in cytRNA-seq and nucRNA- seq. e Percent proportion of abundance of transcripts in the cytoplasm. f Differential expression analysis between cytRNA and nucRNA. Genes enriched in cytRNA are on the right-hand side. Blue, genes with p values less than 0.001 and absolute log2 fold changes greater than unity. g Correlation coefficients of gene expression pattern computed with respect to the conventional scRNA-seq; our novel in silico single-cell normalization showed the best correlation with the scRNA-seq. We also include correlation of nucRNA vs. its in silico single cell than for cytRNA. The yield of cDNA with nucRNA was cDNA from the nucRNA fractions is due to the smaller on par with that of single nuclei prepared with an amount of RNA in a nucleus compared to the cytRNA off-the-shelf kit (PARIS Kit, Thermo Fisher Scientific) in amount for the same cell. The total amount of cDNA which the cell membrane was lysed with a chemical per single cell was 26 ± 16% less than that obtained agent. We thus hypothesize that the smaller amount of with a conventional single-cell protocol on average Abdelmoez et al. Genome Biology (2018) 19:66 Page 4 of 11 (Additional file 1: Figure S2a). We attribute this as with 0.01 < TPM < 1 in nucRNA-seq, bulk nucRNA-seq mainly due to the loss at collecting cytRNA from the detected 4034 genes. This supports the conclusion that outlet well after ITP using a standard micropipette [17]. the detection of genes was not an artifact. We further evaluated this issue by analyzing with the coverages of SINC-seq dissects the difference in subcellular gene the three lowest abundant genes (TPM~ 0.01) expression (Additional file 1:FigureS4b–d). Compared to conventional To benchmark the technical aspects of SINC-seq, we scRNA-seq, the in silico single-cell data showed a correl- assessed the sensitivity and repeatability of gene expres- ation with a mean coefficient of correlation r = 0.866 com- sion analyses with an in silico single-cell analysis. In this puted with a log-transformed expression (log10(TPM + 1)) assessment, we used 56 pairs of nucRNA-seq and (Fig. 1g). The nucRNA-seq also showed a correlation with cytRNA-seq data taken with unperturbed K562 cells scRNA-seq (different cells) with a mean coefficient of correl- which were cultured under standard conditions (without ation r = 0.711, consistent with the results of Grindberg et NaB treatment). (See a comprehensive benchmark of al. [7]. Importantly, nucRNA-seq showed a higher correl- SINC-seq in Additional file 1: Figures S3–S6 and the ation with its in silico single cell (same cells) than those with Supplementary Information section.) SINC-seq consist- single cells (different cells), supporting the conclusion that ently detected 6210 ± 1400 (mean ± SD) and 5690 ± 1500 the nucRNA-seq expression also reflects the transient gene genes per cytRNA and nucRNA, respectively, and 8200 expression of its respective single cell. This finding indicates ± 1100 genes per cell with transcripts per million (TPM) that SINC-seq for the first time reveals the correlation of greater than 1 (see Fig. 1d and comprehensive data in transient gene expression between nucRNA and cytRNA. Additional file 1: Figure S3). SINC-seq also revealed that Combined, an average gene expression profile of 12 in silico ~ 16% of transcripts were in the nucleus and ~ 84% in single cells (see Additional file 1: Figure S7 for the definition the cytoplasm (see Fig. 1e and also section titled “Com- of 12 in silico single cells) showed an excellent matching puting cyt-normalized vs. nuc-normalized data: in silico with that of 12 scRNA-seq (Pearson correlation coefficient single-cell normalization” in Methods) and also showed of r = 0.972, see Additional file 1:FigureS3q).Thetotal enriched expression of 226 and 3035 genes, respectively number of detected genes (TPM > 1) in the 12 in silico (Fig. 1f). Grindberg et al. [7] performed enrichment ana- single cells at an average sequencing depth of 2.2 M reads lyses with individually normalized reads per kilobase per was 11,131, of which 9757 genes were also detected from 12 million mapped reads (RPKM) values from nuclei and scRNA-seq averages (Additional file 1: Figure S3t). We again single cells. Compared to Grindberg’s work, our in silico stress that the in silico normalization and resulting single-cell approach uniquely enables quantitative com- scRNA-seq is a novel method which uniquely leverages a parison of gene expression and enrichment analyses physical separation and recovery as well as minimal between cytoplasm and nucleus. On average, SINC-seq cross-contamination enabled by our electrophoretic displayed about a 5.3% smaller fraction of detected fractionation. genes, 8070 ± 1100 genes at the sum total reads of 2.2 million (M) reads, than the conventional scRNA-seq, which detected 8520 ± 1100 genes (n = 12, 2.2 M reads). Cell-cycle-related genes show correlated expression in Again, we attribute this mainly to the losses associated cytoplasm and nucleus with the collection of the cytRNA from the outlet well To view the landscape of the correlation between nucRNA of the chip using a standard micropipette [17]. Notably, and cytRNA, we computed the cross-correlation of indi- our in silico single-cell data showed a wider dynamic vidual genes as a measure of covariation in the two subcel- range in the detection of genes as compared to lular compartments with 56 pairs of data taken with scRNA-seq (Additional file 1: Figure S4a). We attribute unperturbed K562 cells, and we ranked the genes by order this improvement in the dynamic range of in silico single of the value of the coefficient of correlation (Fig. 2a). We cells to the sequencing depth in nucRNA-seq. SINC-seq identified 1720 positively correlated genes and 29 nega- analyzed RNA expressions in the nucleus and cytoplasm tively correlated genes with p < 0.05. Gene ontology ana- individually with a similar number of sequencing reads. lysis revealed that the correlated genes had cell cycle as an This means that nucRNA-seq has more depth per unit enriched function (Fig. 2b) and the negative correlation RNA template than cytRNA because of the smaller had RNA splicing (Fig. 2c). We note that identification of amount of nuclear RNA. This feature facilitates the uncorrelated genes with this approach was technically detection of low abundant nuclear RNAs, hence result- challenging, as the method resulted in many uncorrelated ing in a larger dynamic range for SINC-seq in silico cell genes with negligible significance (making it impossible to analysis. To assess the validity of the gene detection with neglect the null hypothesis). In this section, we thus TPM < 1 in nucRNA-seq, we evaluated the detection of focused analyses on positively correlated and negatively genes compared with bulk nucRNA-seq. Of 4235 genes correlated genes identifiable with statistical significance. Abdelmoez et al. Genome Biology (2018) 19:66 Page 5 of 11 ab c ef g Fig. 2 Landscape of cross-correlation between cytRNA and nucRNA unveiled transcriptional oscillation of cell-cycle genes in nucRNA highly correlated with expression in cytRNA. a Quantile plot of genes sorted in order of the coefficients of cross-correlation. b Gene ontology analysis with positively correlated genes (p < 0.05) in the quantile plot and c negatively correlated genes. d–f Cell-cycle genes in in silico single-cell data, cytRNA, and nucRNA show correlation with in-phase genes (G1 vs. G1 or G2 vs. G2) and negative correlation with out-of-phase genes (G1 vs. G2). List of genes of G1 and G2 phases are provided in Additional file 1: Figure S8. g Transcriptional oscillation of cell-cycle genes in nucRNA cross- correlated with the gene expression in cytRNA To dissect the correlated gene expression, we focused To explore how the transcriptional oscillations in the on cell cycle based on transcriptional oscillations [2]and nucRNA modulated the gene expression in the cytRNA, we phase-score analysis [3](Additional file 1: Supplementary extended the analyses to compute the cross-correlation Information). Apart from the analysis of the correlation between cytRNA and nucRNA with cell-cycle genes. The landscape, we obtained a list of cell-cycle genes [29]and cell-cycle genes showed synchronized oscillation in each of extended an approach proposed by Klein et al. [2]to the two subcellular compartments (Fig. 2g), consistent with observe the behavior of individual cell-cycle genes. The in the observation that the subpopulations segregated into the silico single-cell data showed the progression of the cell G1 and G2 groups showed corresponding up-regulation cycle with a correlated variation of in-phase genes and and down-regulation of G1 and G2 genes (Additional file 1: negative correlation of out-of-phase genes (G1 vs. G2) Figure S8e, f, and Supplementary Information). Together, (Fig. 2d), consistent with scRNA-seq of K562 [2], and also these results suggest that the cytRNA and nucRNA have with the progression of the phase score (Additional file 1: similar expression patterns of cell-cycle genes and that both Figure S8a). Similarly, both cytRNA-seq and nucRNA-seq of them solely have a potency to detect the cell cycle. data revealed the cell cycle (Fig. 2e, f and Additional file 1: Figure S8b, c). Notably, we found that the negative correl- Nuclear-retained introns attenuate the transcriptional ation among out-of-phase genes was slightly higher in the oscillation − 26 nucRNA (U test, p =8×10 , Additional file 1:Figure We next studied the retained intron (RI)-mediated S8d). On the other hand, the correlation among in-phase regulation of mRNA transport [30–32] leveraging the genes was higher in the cytRNA, which may indicate intron-rich reads with nucRNA-seq of SINC-seq (see further modulation in the cytoplasm. Additional file 1: Figure S9 for comprehensive statistics Abdelmoez et al. Genome Biology (2018) 19:66 Page 6 of 11 on intron detection). After filtering RIs (Methods), genes contained three splicing-related genes [36] and a SINC-seq detected 1760 ± 570 RIs per cell, of which small nucleolar RNA (snoRNA) host gene [37]. These data 1290 ± 540 and 780 ± 290 were detected with lead us to hypothesize that the NRIs likely attenuate the nucRNA-seq and cytRNA-seq, respectively (Fig. 3a). We transcriptional oscillation in the nucleus via fine-tuning of identified 242 nuclear-retained introns (NRIs) in 213 the RNA metabolism in order to maintain the gene genes (Methods). Gene ontology analysis [33] deter- expression and functional RNAs in the cytoplasm. mined that the 213 genes had enriched functions like metabolism of RNA and RNA splicing (Fig. 3b), consist- Sodium butyrate treatment on K562 drives diverging ent with previous studies [31, 34, 35]. We examined the gene expression in subcellular compartments relationship between the probability of NRI and gene To explore the correlation dynamics under perturbation, expression in each individual (subcellular) fraction. We we further performed SINC-seq by differentiating K562 observed a positive correlation between them in the cells along the erythroid lineage with NaB (see nucRNA (r = 0.42, p < 0.01, Fig. 3c), while there was no Methods), sampling 8–13 cells per day over 5 days. We correlation in the cytRNA (r = − 0.006, p = 0.9, Fig. 3d). used NaB to modulate the histone acetylation process In contrast, we observed different enriched functions and observed transient effects of the epigenetic modifi- with cytoplasmic-enriched RIs (CRIs) and no correlation cation on the correlation through the changes in tran- between probabilities of CRI and gene expression scriptional activation and export of mRNA to the (Additional file 1: Figure S10a–c). For a better under- cytoplasm. With 41 successful SINC-seq data sets (82 standing of the function of NRIs, we examined the RNA-seq data in total), we detected differentially expression patterns of the top seven genes that were expressed genes (DEGs); 264 up-regulated and 177 highly associated with NRI-mediated regulation (Fig. 3e–g down-regulated genes were in cytRNA, and 64 and Additional file 1: Figure S10d–h). Notably, the seven up-regulated and 2 down-regulated genes were in a bc ef g Fig. 3 NRI-mediated attenuation of transcriptional oscillation in nucRNA. a Heatmap of the probability of RI in cytRNA and nucRNA fractions. NRI was identified in the upper right region indicated with the broken white line. b Gene ontology analysis with NRI, showing enriched functions like metabolism of RNA and RNA splicing. c, d Correlation analysis between the probability of NRI and the fold change of gene expression among cells with NRI and without NRI (Spl spliced) in nucRNA (upper panel) and cytRNA (lower panel), respectively. e Expression of seven genes that were highly regulated by NRI in nucRNA (p < 0.0003, Mann-Whitney U test), comparing with NRI vs. without NRI (Spl) in an individual fraction. f Coverage of SRSF5 and g HNRNPDL genes showing higher intron reads in the nuclear fraction. Coverages of ARGLU1, GAS5 (SNHG2), FBXO9, VAMP2, and PI4KAP1 genes are provided in Additional file 1: Figure S10 Abdelmoez et al. Genome Biology (2018) 19:66 Page 7 of 11 nucRNA (Additional file 1: Figure S11). We examined cytRNA and nucRNA. To leverage the different behav- the dynamics of the cross-correlation of DEG expression iors of nucRNA and cytRNA, we introduced a between cytRNA and nucRNA along the differentiation localization-embedded principal component analysis by ordering cells with the pseudotime computed using (L-PCA) that computed principal components (PCs) Monocle (version 2.4.0) [38, 39] to in silico single-cell with the subcellular gene expressions of DEGs (Fig. 4b). data (Fig. 4a, Additional file 1: Figure S12a). To provide In the L-PCA, we treated genes in cytRNA and nucRNA the statistical significance, we divided cells into five as different ones and represented a single cell as a vector groups with the pseudotime (see the divisions in Fig. 4a) having double dimension instead of using the conven- and evaluated the cross-correlation of gene expression tional approach. As expected, L-PCA resolved the trajec- between nucRNA and cytRNA (Additional file 1: Figure tory of the differentiation slightly differently compared S12b). The cross-correlation between cytRNA and to a conventional PCA that computed PCs with in silico nucRNA exhibited a gradual decrease with increasing single-cell data (Fig. 4c). To further corroborate the pseudotime, suggesting that the subcellular gene expres- L-PCA, we performed PCA on data sets of an individual sion patterns behaved differently and diverged as the fraction, showing this specific and unique transient in differentiation proceeded. We also found that non-DEGs the nucRNA — that is the day 3 cluster was furthest showed a less significant change in the coefficient of from the day 0 cluster in nucRNA (Additional file 1: Fig- cross-correlation with the differentiation (Additional file 1: ure S12f, g). With our data set, it was difficult to con- Figure S12c). We further explored this observation with clude whether L-PCA can offer better clustering than the control experiments (Additional file 1: Figure S12c–e). conventional PCA. However, these results at least dem- We found that the cross-correlation with DEGs showed onstrated that the L-PCA with SINC-seq is practical and no apparent change with unperturbed cells along the potentially useful when analyzing a biological process in- pseudotime (Additional file 1: Figure S12d). These find- volving regulations at multiple layers of single cells. ings and data with population controls at day 0 and day 4 (Additional file 1: Figure S12e) indicate that the change in Discussion the cross-correlation was not an artifact but reflects the A fundamental question is how the transcriptional oscil- dynamics of cytRNA and nucRNA along the differenti- lation in the nucRNA, which is inherently stochastic, is ation. Notably, the cross-correlation between nucRNA on transported to and correlated with gene expression in day 4 and cytRNA on day 1 (near the top right corner in the cytRNA. SINC-seq enabled direct and quantitative Fig. 4a) was lower compared to that between nucRNA on comparison of gene expressions between a nucleus, a day 1 and cytRNA on day 4 (near the bottom left corner), cytoplasm, and a whole cell of the same single cell. This suggesting that nucRNA was the driver of the diverging comparison reveals that the cell may conceivably gene expression. fine-tune a portion of its expression upon transport to We hypothesize that the divergence of gene expres- the cytoplasm (e.g., with NRI genes), while preserving sions in the two subcellular compartments reflects the correlation of other portions of its expression upon transient responses of different regulatory pathways in transport (e.g., with cell-cycle-related genes). SINC-seq nuc Pseudotime bc 01 2 3 4 L-PCA Conventional PCA (51.0) 0.8 Day4 Day2 Day1 (32.5) 10 (66.8) Day3 0.6 Day0 Day2 Day0 0.4 Day1 −10 Day3 Day4 (67.4) −20 0.2 −20 −30 −25 0 25 −25 0 25 PC1 PC1 Fig. 4 Differentiation of K562 cells to erythroid cells shows a dynamical change of cross-correlation between cytRNA and nucRNA. a The cross- correlation of DEG expression between cytRNA and nucRNA along with the pseudotime. Numbers along axes show groups of cells used in the analysis shown in Additional file 1: Figure S12b. b L-PCA analysis of differentiating K562 cells to erythroid cells by sodium butyrate treatment. c Conventional PCA analysis with differentiating K562 cells cyt Pseudotime 43 2 1 0 PC2 PC2 Abdelmoez et al. Genome Biology (2018) 19:66 Page 8 of 11 also revealed that the cells under external perturbation Buffers dynamically alter the correlation and exhibit a unique We designed buffers for ITP-based selective extraction, trajectory of differentiation at subcellular resolution. separation (from the trapped nucleus), purification, and These findings shed new light on the characteristics of transport of cytRNA to the cytRNA output well of the post-transcriptional regulations with a single cell and chip (see more details in Shintaku et al. [16] and subcellular compartment resolution. Kuriyama et al. [17]). The leading electrolyte buffer (LE) Our study also suggests a compelling caution to an components were 50 mM Tris and 25 mM HCl contain- approach that approximates the transcriptomic profile of ing 0.4% poly(vinylpyrrolidone) (PVP) (calculated pH of the whole cell with that of a single compartment without 8.1). The trailing electrolyte buffer (TE) components validation. The SINC-seq platform will be broadly were 50 mM imidazole and 25 mM HEPES containing applicable to different types of cells as long as they are (initial calculated pH of 8.3) 0.4% PVP. We included isolated as singles. The method thus will contribute to PVP to suppress electroosmotic flow. We purchased validate existing subcellular RNA-seq methods [4, 5, 9– Tris, HEPES, imidazole, and HCl from Sigma-Aldrich, 11, 19, 20] and also define their limitations. and PVP (molecular weight 1 MDa) from Polyscience. We prepared all solutions in UltraPure DNase-/RNase-free deionized (DI) water (Life Technologies). Conclusion To dissect transcriptional correlation in subcellular com- Microfluidic system setup partments, we devised SINC-seq, which enables inte- We fabricated polydimethylsiloxane (PDMS, Sylgard grated nuclear and cytoplasmic RNA-seq of single cells 184, Dow Corning) microchannel superstructures by coupling physical fractionation of cytRNA from the (Additional file 1: Figure S1a, b) with a soft lithography nucleus of a single cell with a high-throughput and bonded them to a glass substrate. SU-8 (SU-8 2025, RNA-seq. Leveraging SINC-seq, we explored the land- MicroChem) molds were prepared on glass substrates scape of the correlation between nucRNA and cytRNA with the microchannel patterns made of chromium thin with a total of 84 K562 cells, which corresponds to 168 films, exposing the SU-8 to ultraviolet light through the RNA-seq libraries. The SINC-seq data unveiled three pattern. The nominal channel width and depth of the distinct natures of correlation among cytRNA and microchannels were 50 μm and 25 μm, respectively. We nucRNA that reflected the physiological state of single designed 3-μm-wide and 5-μm-long hydrodynamic traps. cells: highly correlated expression in cell-cycle-related We have optimized the size of the hydrodynamic trap as genes, the distorted correlation via NRIs, and the correl- 3 μm wide and 25 μm deep so that it can capture a ation dynamics along the differentiation of K562 cells to single cell and trap a nucleus during the electric field erythroid cells under sodium butyrate perturbation. extraction of cytRNA. We observed some deformation These data provide unique insights into the regulatory of cells due to the presence of the pressure-driven flow, network of mRNA from the nucleus toward the cyto- but observed no mechanical lyses prior to the applica- plasm at the single-cell level. tion of the electric field. This trapping process had a timescale of several minutes. Methods Before each experiment, we preconditioned the micro- Cells channel by filling the inlet and outlet wells with washing We purchased K562 cells (human lymphoblast, solutions and applying a vacuum to the waste well. Our chronic myelogenous leukemia) from RIKEN BioRe- washing process was as follows: 1 M NaOH for 1 min, source Center and the Japanese Collection of 1 M HCl for 1 min, and DI water for 1 min. All washing Research Bioresources (JCRB) cell bank. We cultured solutions contained 0.1% Triton X-100 to suppress the K562 cells in RPMI-1640 Medium (Life Technolo- bubble clogging in the hydrophobic microchannel. gies) with 10% fetal bovine serum and 1% Following this, we loaded 9.5 μL of LE and TE to the penicillin-streptomycin-glutamine at 37 °C in 5% CO . outlet and inlet wells, respectively, and briefly applied a We washed the cells with phosphate-buffered saline vacuum to the waste well to exchange the solution in once and suspended them in a sample buffer the microchannel with LE and TE. The hydrodynamic containing 50 mM imidazole, 25 mM 4-(2-hydro- pressure induced by buffers in the inlet and outlet wells xyethyl)-1-piperazineethanesulfonic acid (HEPES), and created a pressure-driven laminar flow from both inlet 175 mM sucrose (pH 7.6) at a concentration of ~ 0.8 and outlet wells toward the waste well and formed a cells/μL and stored them on ice until the experiments stable LE-TE interface at the junction of three micro- were performed. To differentiate the K562 cells, we channels. We then picked up a single cell of interest incubated them with 1 mM NaB (Sigma-Aldrich, from the cell suspension by aspirating 1 μL using a B5887) and harvested them after 96 h of induction. standard micropipette. We observed this process using a Abdelmoez et al. Genome Biology (2018) 19:66 Page 9 of 11 microscope and confirmed aspiration of a single cell. using a Nextera XT DNA Sample Preparation Kit (Illu- We then dispensed this 1-μL volume into the inlet well mina) and purified the cDNA library following the man- and introduced it into the microchannel via the ufacturer’s protocol, except that we eluted the cDNA pressure-driven flow. Once we visually confirmed the sample with 24 μL instead of 50 μL (see Additional file 1: captured single cell at the hydrodynamic trap (Fig. 1c), Figure S2a for yields of cDNA). We pooled 98–108 we added 9.5 μL of the TE to the waste well to reduce libraries and sequenced them on an Illumina HiSeq the pressure-driven flow. We aborted our protocol in 2500 with 100-base paired-end reads to an average depth cases where we observed two or more cells at the hydro- of 4.64 million reads (Additional file 1: Figure S2b, c). dynamic trap. We placed 300-μm-diameter platinum We mapped the trimmed sequencing reads to the tran- wire electrodes into the wells and applied − 150 V, − scripts derived from the human reference genome 170 V, and 0 V to the electrodes at the inlet, waste, and (GRCh37.75) using the STAR (version 2.5.1b) mapping outlet wells, respectively. The DC voltage created a con- program [40] with ENCODE options, and calculated centrated electric field at the hydrodynamic trap expression estimates with TPM using RNA-Seq by (Additional file 1: Figure S1d) and lysed the cytoplasmic Expectation-Maximization (RSEM v1.3.0) [41]. membrane within 1 s. Appropriate placement of the ITP We performed DE analyses individually with cytRNA buffers with the DC electric field enabled an immediate and nucRNA comparing day 1–day 4 samples to day 0 transition from the lysis to an ITP process that collects samples. We used the MATLAB functions “nbintest” and focuses cytRNA into an ITP-zone, TE-to-LE inter- and “mafdr” and identified significance with p values less face. At 40 s, we changed the voltages to − 350 V and − than 0.001 and absolute log2 fold changes greater than 510 V at the inlet and waste wells, respectively, to accel- unity. erate the migration of the ITP zone. The ITP zone trans- ported the cytRNA to the output well in about 100 s, Analysis of intron retention while the nucleus was retained at the hydrodynamic We computed intron expressions with fragments per trap. We also monitored the current during the extrac- kilobase of intron per million mapped reads (FPKM) tion with a computer running a custom MATLAB using 347,041 unique introns (longer than 50 nt) on the (Mathworks, Inc., Natick, MA, USA) script. The magni- genome annotation with 56 SINC-seq data of K562 cells tude of the current decreased as the ITP zone (contain- under standard culturing conditions. On average, ing the focused cytRNA) advanced in the channel and as SINC-seq yielded 14.8% reads mapped to introns with the lower conductivity TE replaced the higher conduct- nucRNA-seq, but only 1.1% with cytRNA-seq. Further, ivity LE (Additional file 1: Figure S1c). The current SINC-seq detected 38,800 ± 10,000 and 34,800 ± 8000 signal plateaued near t = 100 s, coincident with the time unique introns in nucRNA-seq and cytRNA-seq, re- at which the focused cytRNA eluted into the outlet well. spectively, and 56,300 ± 9500 per cell with FPKM of We deactivated the voltages at 200 s and used a standard more than 0. We identified an RI that had at least 10% pipette to transfer two aliquots from the chip: 9.5 μL expression level of the gene, 95% coverage in the in- from the outlet well containing the cytRNA and 1 μL tronic region, and non-zero expression in the adjacent containing the cell nucleus from the inlet well. Detailed exon. We discarded intron reads locating on a gene with descriptions of a similar protocol and chip, together with less than 2 TPM. On the other hand, we identified a a narrated video description, were reported by Kuriyama fully spliced intron that had less than 1% expression of et al. [18]. the gene and 50% coverage in the adjacent exon. We dis- carded intron reads that failed the preceding criteria. Library preparation and mapping analysis We validated the RI identification with splice site scores We synthesized respective cDNA libraries from the [42], which showed lower values with RIs than fully fractionated cytRNA and nucRNA separately using spliced introns, using 9mer (exonic 3mer + intronic − 16 Smart-seq2 (SMART-seq v4 Ultra Low Input RNA Kit 6mer) around the 5′ splice site (p value < 2.2 × 10 , U for Sequencing, Clontech) with 18 polymerase chain test), and 23mer (intronic 20mer + exonic 3mer) around − 12 reaction (PCR) cycles followed by purification with the 3′ splice site (p value < 2.5 × 10 , U test). In total, Agencourt AMPure XP (Beckman Coulter). We exam- we detected 17,277 RIs, of which 14,134 and 2316 RIs ined the yield and quality of cDNA, respectively, with a had a higher probability of RI in nucRNA and cytRNA, Qubit 2.0 Fluorometer (Thermo Fisher Scientific) and respectively. The RI enrichment in the nucleus was quantitative PCR (qPCR) targeting GAPDH (glyceralde- consistent with that for previous studies [30–32]. hyde-3-phosphate dehydrogenase, Hs02758991_g1, To identify NRIs, we calculated the probability of in- Thermo Fisher Scientific) and HBG (gamma-globin tron retention defined as the proportion of cells with the genes, Hs00361131_g1, Thermo Fisher Scientific). We RI and identified NRIs that had 0.25 higher probability performed the tagmentation reaction with 200 pg cDNA in the nuclear fraction than in the cytoplasmic fraction. Abdelmoez et al. Genome Biology (2018) 19:66 Page 10 of 11 We then filtered unique NRIs, discarding smaller NRIs Additional files that had overlap with a long NRI. On the other hand, we Additional file 1: Supplementary information and Figures S1–S12. identified CRIs that had 0.25 higher probability in the (DOCX 4587 kb) cytoplasmic fraction than in the nuclear fraction. Additional file 2: Movie S1. Electrical lysis and RNA extraction visualized by SYBR Green II. (MOV 1279 kb) Computing cyt-normalized vs. nuc-normalized data: in Additional file 3: Table S1. Quality control of SINC-seq samples. silico single-cell normalization (XLSX 21 kb) We computed in silico single-cell RNA-seq (scRNA-seq) Additional file 4: Movie S2. Experimental run #55 with unsuccessful fractionation. (MP4 1372 kb) data with cytRNA-seq and nucRNA-seq data, scaling the Additional file 5: Movie S3. Experimental run #69 with unsuccessful raw TPM values and combining the cytRNA-seq and fractionation. (MP4 1116 kb) nucRNA-seq data as Acknowledgements TPM ¼ TPM þ TPM ; ðS1Þ in−silico cyt nuc We acknowledge Yuka Ozaki of Kyoto University for technical support in the preparation of the sequencing libraries. TPM ¼ αTPM ; ðS2Þ cyt Funding This work was funded by the ImPACT Program of the Council for Science, TPM ¼ βTPM ; ðS3Þ nuc Technology, and Innovation (Cabinet Office, Government of Japan) and also by the Japan Society for the Promotion of Science under 26289035 and where α and β are normalization factors, which make the summation of the TPM values, TPM , of the in in silico Availability of data and materials silico single-cell data to be one million. We here write The data from this study have been deposited into the Sequencing Read Archive and are available for download under accession SRP119800 [43]. raw TPM values with an asterisk and TPM values of in silico single-cell data, cytRNA-seq, and nucRNA-seq, Authors’ contributions respectively, with their subscripts. We calculate α and β MNA performed the experiments; YO, KI, HS, and HN analyzed the data. SU, RY, and HK helped to analyze the data. HS and JGS designed and supervised using ΔCt of qPCR data taken at the QC of cDNA the project. All authors wrote the manuscript. All authors read and approved (Additional file 3: Table S1) as the final manuscript. ΔCt 2 Ethics approval α ¼ ; ðS4Þ Does not apply to this study. TPM cyt;geneA ΔCt þ 2 TPM nuc;geneA Competing interests H.S., Y.O., and S.U. submitted patent applications (Japanese Patent Application No. 2016-014708, No. 2016-177,163, and PCT/JP2017/003069) β ¼ 1−α; ðS5Þ relating to this work. Correspondence and requests for materials should be addressed to H.S. ([email protected]). * * where TPM , TPM , and ΔCt are, respect- cyt,geneA nuc, geneA ively, the raw TPM value of gene A with cytRNA-seq, Publisher’sNote Springer Nature remains neutral with regard to jurisdictional claims in the raw TPM value of gene A with nucRNA-seq, and published maps and institutional affiliations. ΔCt = Ct – Ct which is ΔCt with nuc, geneA cyt, geneA, respect to gene A. We calculated pairs of α and β with Author details Department of Micro Engineering, Graduate School of Engineering, Kyoto GAPDH and HBG genes (HBG1 + HBG2)as gene As, University, Kyoto, Japan. Microfluidics RIKEN Hakubi Research Team, RIKEN and used the mean α and β to compute the TPM. Cluster for Pioneering Research, Saitama, Japan. Medical Research Support We here note that α and β, respectively, indicate the Center, Graduate School of Medicine, Kyoto University, Kyoto, Japan. Department of Biological Sciences, Graduate School of Science, The (complementary) fractions of cytoplasmic transcripts University of Tokyo, Tokyo, Japan. Department of Hematology, Faculty of and nuclear transcripts in the in silico single cells. We Medicine, University of Tsukuba, Tsukuba, Japan. Department of Mechanical thus can quantify the fraction of the cytoplasmic tran- Engineering, Stanford University, Stanford, CA, USA. script as shown in Fig. 1e. Received: 13 October 2017 Accepted: 7 May 2018 L-PCA The PCA for the conventional scRNA-seq uses an (n × m) References 1. Wu AR, et al. Quantitative assessment of single-cell RNA-sequencing matrix of gene expressions (n = number of genes) with methods. Nat Methods. 2014;11(1):41–6. multiple samples (m = number of samples); however, our 2. Klein AM, et al. Droplet barcoding for single-cell transcriptomics applied to L-PCA uses a (2n × m)matrix of geneexpressions,having embryonic stem cells. Cell. 2015;161(5):1187–201. 3. Macosko EZ, et al. Highly parallel genome-wide expression profiling of a twofold dimension compared to the conventional individual cells using nanoliter droplets. Cell. 2015;161(5):1202–14. matrix, derived from cytRNA-seq and nucRNA-seq. The 4. Habib N, et al. Div-Seq: single-nucleus RNA-Seq reveals dynamics of rare L-PCA was performed using PCA with “prcomp” in R. adult newborn neurons. Science. 2016;353(6302):925–8. Abdelmoez et al. Genome Biology (2018) 19:66 Page 11 of 11 5. Lake BB, et al. Neuronal subtypes and diversity revealed by single-nucleus 37. Williams GT, Farzaneh F. Are snoRNAs and snoRNA host genes new players RNA sequencing of the human brain. Science. 2016;352(6293):1586–90. in cancer? Nat Rev Cancer. 2012;12(2):84–8. 6. Tanay A, Regev A. Scaling single-cell genomics from phenomenology to 38. Qiu X, et al. Single-cell mRNA quantification and differential analysis with mechanism. Nature. 2017;541(7637):331–8. Census. Nat Methods. 2017;14(3):309–15. 7. Grindberg RV, et al. RNA-sequencing from single nuclei. Proc Natl Acad Sci 39. Trapnell C, et al. The dynamics and regulators of cell fate decisions are U S A. 2013;110(49):19802–7. revealed by pseudotemporal ordering of single cells. Nat Biotechnol. 2014; 32(4):381–6. 8. Krishnaswami SR, et al. Using single nuclei for RNA-seq to capture the 40. Dobin A, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. transcriptome of postmortem neurons. Nat Protoc. 2016;11(3):499–524. 2013;29(1):15–21. 9. Lacar B, et al. Nuclear RNA-seq of single neurons reveals molecular 41. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data signatures of activation. Nat Commun. 2016;7:11022. with or without a reference genome. BMC Bioinformatics. 2011;12:323. 10. Habib N, et al. Massively parallel single-nucleus RNA-seq with DroNc-seq. 42. Yeo G, Burge CB. Maximum entropy modeling of short sequence motifs with Nat Methods. 2017;14(10):955–8. applications to RNA splicing signals. J Comput Biol. 2004;11(2-3):377–94. 11. Gao R, et al. Nanogrid single-nucleus RNA sequencing reveals phenotypic 43. Abdelmoez MN, et al., Gene expressions in nucleus and cytoplasm at the diversity in breast cancer. Nat Commun. 2017;8(1):228. single cell resolution. Sequencing Read Archive. 2017. https://www.ncbi.nlm. 12. Macaulay IC, et al. G&T-seq: parallel sequencing of single-cell genomes and nih.gov/bioproject/PRJNA413830. transcriptomes. Nat Methods. 2015;12(6):519–22. 13. Angermueller C, et al. Parallel single-cell sequencing links transcriptional and epigenetic heterogeneity. Nat Methods. 2016;13(3):229–32. 14. Cheow LF, et al. Single-cell multimodal profiling reveals cellular epigenetic heterogeneity. Nat Methods. 2016;13(10):833–6. 15. Han L, et al. Co-detection and sequencing of genes and transcripts from the same single cells facilitated by a microfluidics platform. Sci Rep. 2014;4:6485. 16. Shintaku H, et al. On-chip separation and analysis of RNA and DNA from single cells. Anal Chem. 2014;86(4):1953–7. 17. Kuriyama K, Shintaku H, Santiago JG. Isotachophoresis for fractionation and recovery of cytoplasmic RNA and nucleus from single cells. Electrophoresis. 2015;36(14):1658–62. 18. Kuriyama K, Shintaku H, Santiago JG. Protocol for microfluidic system to automate the preparation and fractionation of the nucleic acids in the cytoplasm versus nuclei of single cells. Bio-protocol. 2016;6(12):e1844. 19. Hou Y, et al. Single-cell triple omics sequencing reveals genetic, epigenetic, and transcriptomic heterogeneity in hepatocellular carcinomas. Cell Res. 2016;26(3):304–19. 20. Hu Y, et al. Simultaneous profiling of transcriptome and DNA methylome from a single cell. Genome Biol. 2016;17:88. 21. Kohler A, Hurt E. Exporting RNA from the nucleus to the cytoplasm. Nat Rev Mol Cell Biol. 2007;8(10):761–73. 22. Gerstberger S, Hafner M, Tuschl T. A census of human RNA-binding proteins. Nat Rev Genet. 2014;15(12):829–45. 23. Wang J, Lu C, et al. Kinetics of NF-κB nucleocytoplasmic transport probed by single-cell screening without imaging. Lab Chip. 2010;10(21):2911–6. 24. Zhan Y, Lu C, et al. One-step extraction of subcellular proteins from eukaryotic cells. Lab Chip. 2010;10(16):2046–8. 25. Qu, Marshall LA, Santiago JG. Simultaneous purification and fractionation of nucleic acids and proteins from complex samples using bidirectional isotachophoresis. Anal Chem. 2014;86(15):7264–8. 26. Picelli S, et al. Smart-seq2 for sensitive full-length transcriptome profiling in single cells. Nat Methods. 2013;10(11):1096–8. 27. Ziegenhain C, et al. Comparative analysis of single-cell RNA sequencing methods. Mol Cell. 2017;65(4):631–43. e4 28. Habib N, et al. DroNc-Seq: deciphering cell types in human archived brain tissues by massively-parallel single nucleus RNA-seq. bioRxiv. 2017; https:// doi.org/10.1101/115196 29. Whitfield ML, et al. Identification of genes periodically expressed in the human cell cycle and their expression in tumors. Mol Biol Cell. 2002;13(6):1977–2000. 30. Ninomiya K, Kataoka N, Hagiwara M. Stress-responsive maturation of Clk1/4 pre-mRNAs promotes phosphorylation of SR splicing factor. J Cell Biol. 2011;195(1):27–40. 31. Boutz PL, Bhutkar A, Sharp PA. Detained introns are a novel, widespread class of post-transcriptionally spliced introns. Genes Dev. 2015;29(1):63–80. 32. Naro C, et al. An orchestrated intron retention program in meiosis controls timely usage of transcripts during germ cell differentiation. Dev Cell. 2017; 41(1):82–93. e4 33. Tripathi S, et al. Meta- and orthogonal integration of influenza “OMICs” data defines a role for UBR4 in virus budding. Cell Host Microbe. 2015;18(6):723–35. 34. Iida K, Go M. Survey of conserved alternative splicing events of mRNAs encoding SR proteins in land plants. Mol Biol Evol. 2006;23(5):1085–94. 35. Hsu TY, et al. The spliceosome is a therapeutic vulnerability in MYC-driven cancer. Nature. 2015;525(7569):384–8. 36. Lareau LF, et al. Unproductive splicing of SR genes associated with highly conserved and ultraconserved DNA elements. Nature. 2007;446(7138):926–9.

Journal

Genome Biology – Springer Journals

Published: Jun 6, 2018

There are no references for this article.