Single-cell RNA sequencing of a European and an African lymphoblastoid cell line

Private information leakage from single-cell count matrices

The increase in publicly available human single-cell datasets, encompassing millions of cells from many donors, has significantly enhanced our understanding of complex biological processes. However, the accessibility of these datasets raises significant privacy concerns. Due to the inherent noise in single-cell measurements and the scarcity of population-scale single-cell datasets, recent private information quantification studies have focused on bulk gene expression data sharing. To address this gap, we demonstrate that individuals in single-cell gene expression datasets are vulnerable to linking attacks, where attackers can infer their sensitive phenotypic information using publicly available tissue or cell-type-specific expression quantitative trait loci (eQTLs) information. We further develop a method for genotype prediction and genotype-phenotype linking that remains effective without relying on eQTL information. We show that variants from one study can be exploited to uncover private information about individuals in another study.

Evaluating the role of the nuclear microenvironment in gene function by population-based modeling

The nuclear folding of chromosomes relative to nuclear bodies is an integral part of gene function. Here, we demonstrate that population-based modeling-from ensemble Hi-C data-provides a detailed description of the nuclear microenvironment of genes and its role in gene function. We define the microenvironment by the subnuclear positions of genomic regions with respect to nuclear bodies, local chromatin compaction, and preferences in chromatin compartmentalization. These structural descriptors are determined in single-cell models, thereby revealing the structural variability between cells. We demonstrate that the microenvironment of a genomic region is linked to its functional potential in gene transcription, replication, and chromatin compartmentalization. Some chromatin regions feature a strong preference for a single microenvironment, due to association with specific nuclear bodies in most cells. Other chromatin shows high structural variability, which is a strong indicator of functional heterogeneity. Moreover, we identify specialized nuclear microenvironments, which distinguish chromatin in different functional states and reveal a key role of nuclear speckles in chromosome organization. We demonstrate that our method produces highly predictive three-dimensional genome structures, which accurately reproduce data from a variety of orthogonal experiments, thus considerably expanding the range of Hi-C data analysis.

Conformational analysis of chromosome structures reveals vital role of chromosome morphology in gene function

The 3D conformations of chromosomes are highly variant and stochastic between single cells. Recent progress in multiplexed 3D FISH imaging, single cell Hi-C and genome structure modeling allows a closer analysis of the structural variations of chromosomes between cells to infer the functional implications of structural heterogeneity. Here, we introduce a two-step dimensionality reduction method to classify a population of single cell 3D chromosome structures, either from simulation or imaging experiment, into dominant conformational clusters with distinct chromosome morphologies. We found that almost half of all structures for each chromosome can be described by 5-10 dominant chromosome morphologies, which play a fundamental role in establishing conformational variation of chromosomes. These morphologies are conserved in different cell types, but vary in their relative proportion of structures. Chromosome morphologies are distinguished by the presence or absence of characteristic chromosome territory domains, which expose some chromosomal regions to varying nuclear environments in different morphologies, such as nuclear positions and associations to nuclear speckles, lamina, and nucleoli. These observations point to distinct functional variations for the same chromosomal region in different chromosome morphologies. We validated chromosome conformational clusters and their associated subnuclear locations with data from DNA-MERFISH imaging and single cell sci-HiC data. Our method provides an important approach to assess the variation of chromosome structures between cells and link differences in conformational states with distinct gene functions.

A comprehensive single cell data analysis of lymphoblastoid cells reveals the role of super-enhancers in maintaining EBV latency

We probed the lifecycle of Epstein-Barr virus (EBV) on a cell-by-cell basis using single cell RNA sequencing (scRNA-seq) data from nine publicly available lymphoblastoid cell lines (LCLs). While the majority of LCLs comprised cells containing EBV in the latent phase, two other clusters of cells were clearly evident and were distinguished by distinct expression of host and viral genes. Notably, both were high expressors of EBV LMP1/BNLF2 and BZLF1 compared to another cluster that expressed neither gene. The two novel clusters differed from each other in their expression of EBV lytic genes, including glycoprotein gene GP350. The first cluster, comprising GP350- LMP1hi cells, expressed high levels of HIF1A and was transcriptionally regulated by HIF1-α. Treatment of LCLs with Pevonedistat, a drug that enhances HIF1-α signaling, markedly induced this cluster. The second cluster, containing GP350+ LMP1hi cells, expressed EBV lytic genes. Host genes that are controlled by super-enhancers (SEs), such as transcription factors MYC and IRF4, had the lowest expression in this cluster. Functionally, the expression of genes regulated by MYC and IRF4 in GP350+ LMP1hi cells were lower compared to other cells. Indeed, induction of EBV lytic reactivation in EBV+ AKATA reduced the expression of these SE-regulated genes. Furthermore, CRISPR-mediated perturbation of the MYC or IRF4 SEs in LCLs induced the lytic EBV gene expression, suggesting that host SEs and/or SE target genes are required for maintenance of EBV latency. Collectively, our study revealed EBV-associated heterogeneity among LCLs that may have functional consequence on host and viral biology.

Promoter sequence and architecture determine expression variability and confer robustness to genetic variants

Genetic and environmental exposures cause variability in gene expression. Although most genes are affected in a population, their effect sizes vary greatly, indicating the existence of regulatory mechanisms that could amplify or attenuate expression variability. Here, we investigate the relationship between the sequence and transcription start site architectures of promoters and their expression variability across human individuals. We find that expression variability can be largely explained by a promoter's DNA sequence and its binding sites for specific transcription factors. We show that promoter expression variability reflects the biological process of a gene, demonstrating a selective trade-off between stability for metabolic genes and plasticity for responsive genes and those involved in signaling. Promoters with a rigid transcription start site architecture are more prone to have variable expression and to be associated with genetic variants with large effect sizes, while a flexible usage of transcription start sites within a promoter attenuates expression variability and limits genotypic effects. Our work provides insights into the variable nature of responsive genes and reveals a novel mechanism for supplying transcriptional and mutational robustness to essential genes through multiple transcription start site regions within a promoter.

Epstein-Barr virus perpetuates B cell germinal center dynamics and generation of autoimmune-associated phenotypes in vitro

Human B cells encompass functionally diverse lineages and phenotypic states that contribute to protective as well as pathogenic responses. Epstein-Barr virus (EBV) provides a unique lens for studying heterogeneous B cell responses, given its adaptation to manipulate intrinsic cell programming. EBV promotes the activation, proliferation, and eventual outgrowth of host B cells as immortalized lymphoblastoid cell lines (LCLs) in vitro, which provide a foundational model of viral latency and lymphomagenesis. Although cellular responses and outcomes of infection can vary significantly within populations, investigations that capture genome-wide perspectives of this variation at single-cell resolution are in nascent stages. We have recently used single-cell approaches to identify EBV-mediated B cell heterogeneity in de novo infection and within LCLs, underscoring the dynamic and complex qualities of latent infection rather than a singular, static infection state. Here, we expand upon these findings with functional characterizations of EBV-induced dynamic phenotypes that mimic B cell immune responses. We found that distinct subpopulations isolated from LCLs could completely reconstitute the full phenotypic spectrum of their parental lines. In conjunction with conserved patterns of cell state diversity identified within scRNA-seq data, these data support a model in which EBV continuously drives recurrent B cell entry, progression through, and egress from the Germinal Center (GC) reaction. This "perpetual GC" also generates tangent cell fate trajectories including terminal plasmablast differentiation, which constitutes a replicative cul-de-sac for EBV from which lytic reactivation provides escape. Furthermore, we found that both established EBV latency and de novo infection support the development of cells with features of atypical memory B cells, which have been broadly associated with autoimmune disorders. Treatment of LCLs with TLR7 agonist or IL-21 was sufficient to generate an increased frequency of IgD-/CD27-/CD23-/CD38+/CD138+ plasmablasts. Separately, de novo EBV infection led to the development of CXCR3+/CD11c+/FCRL4+ B cells within days, providing evidence for possible T cell-independent origins of a recently described EBV-associated neuroinvasive CXCR3+ B cell subset in patients with multiple sclerosis. Collectively, this work reveals unexpected virus-driven complexity across infected cell populations and highlights potential roles of EBV in mediating or priming foundational aspects of virus-associated immune cell dysfunction in disease.

Time-resolved transcriptomes reveal diverse B cell fate trajectories in the early response to Epstein-Barr virus infection

Epstein-Barr virus infection of B lymphocytes elicits diverse host responses via well-adapted transcriptional control dynamics. Consequently, this host-pathogen interaction provides a powerful system to explore fundamental processes leading to consensus fate decisions. Here, we use single-cell transcriptomics to construct a genome-wide multistate model of B cell fates upon EBV infection. Additional single-cell data from human tonsils reveal correspondence of model states to analogous in vivo phenotypes within secondary lymphoid tissue, including an EBV+ analog of multipotent activated precursors that can yield early memory B cells. These resources yield exquisitely detailed perspectives of the transforming cellular landscape during an oncogenic viral infection that simulates antigen-induced B cell activation and differentiation. Thus, they support investigations of state-specific EBV-host dynamics, effector B cell fates, and lymphomagenesis. To demonstrate this potential, we identify EBV infection dynamics in FCRL4+/TBX21+ atypical memory B cells that are pathogenically associated with numerous immune disorders.

Uncovering the Principles of Genome Folding by 3D Chromatin Modeling

Our understanding of how genomic DNA is tightly packed inside the nucleus, yet is still accessible for vital cellular processes, has grown dramatically over recent years with advances in microscopy and genomics technologies. Computational methods have played a pivotal role in the structural interpretation of experimental data, which helped unravel some organizational principles of genome folding. Here, we give an overview of current computational efforts in mechanistic and data-driven 3D chromatin structure modeling. We discuss strengths and limitations of different methods and evaluate the added value and benefits of computational approaches to infer the 3D structural and dynamic properties of the genome and its underlying mechanisms at different scales and resolution, ranging from the dynamic formation of chromatin loops and topological associated domains to nuclear compartmentalization of chromatin and nuclear bodies.

Discovery and Functional Characterization of Two Regulatory Variants Underlying Lupus Susceptibility at 2p13.1

Genome-wide association studies have identified 2p13.1 as a prominent susceptibility locus for systemic lupus erythematosus (SLE)—a complex, multisystem autoimmune disease. However, the identity of underlying causal variant (s) and molecular mechanisms for increasing disease susceptibility are poorly understood. Using meta-analysis (cases = 10,252, controls = 21,604) followed by conditional analysis, bioinformatic annotation, and eQTL and 3D-chromatin interaction analyses, we computationally prioritized potential functional variants and subsequently experimentally validated their effects. Ethnicity-specific meta-analysis revealed striking allele frequency differences between Asian and European ancestries, but with similar odds ratios. We identified 20 genome-wide significant (p < 5 × 10−8) variants, and conditional analysis pinpointed two potential functional variants, rs6705628 and rs2272165, likely to explain the association. The two SNPs are near DGUOK, mitochondrial deoxyguanosine kinase, and its associated antisense RNA DGUOK-AS1. Using luciferase reporter gene assays, we found significant cell type- and allele-specific promoter activity at rs6705628 and enhancer activity at rs2272165. This is supported by ChIP-qPCR showing allele-specific binding with three histone marks (H3K27ac, H3K4me3, and H3K4me1), RNA polymerase II (Pol II), transcriptional coactivator p300, CCCTC-binding factor (CTCF), and transcription factor ARID3A. Transcriptome data across 28 immune cell types from Asians showed both SNPs are cell-type-specific but only in B-cells. Splicing QTLs showed strong regulation of DGUOK-AS1. Genotype-specific DGOUK protein levels are supported by Western blots. Promoter capture Hi-C data revealed long-range chromatin interactions between rs2272165 and several nearby promoters, including DGUOK. Taken together, we provide mechanistic insights into how two noncoding variants underlie SLE risk at the 2p13.1 locus.

Transcriptome Analysis on Key Metabolic Pathways in Rhodotorula mucilaginosa Under Pb(II) Stress

Rhodotorula mucilaginosa shows adaption to a broad range of Pb²⁺ stress. In this study, three key pathways, i.e., glycolysis (EMP), the tricarboxylic acid (TCA) cycle, and oxidative phosphorylation (OXPHOS), were investigated under 0-2,500 mg · L^-1 Pb stress, primarily based on biochemical analysis and RNA sequencing. R. mucilaginosa cells showed similar metabolic response to low/medium (500/1000 mg · L^-1) Pb²⁺ stress. High (2,500 mg · L^-1) Pb²⁺ stress exerted severe cytotoxicity to R. mucilaginosa. The downregulation of HK under low-medium Pb²⁺ suggested a correlation with the low hexokinase enzymatic activity in vivo. However, IDH3, regulating a key step of circulation in TCA, was upregulated to promote ATP feedstock for downstream OXPHOS. Then, through activation of complex I & IV in the electron transport chain (ETC) and ATP synthase, ATP production was finally enhanced. This mechanism enabled fungal cells to compensate for ATP consumption under low-medium Pb²⁺ toxicity. Hence, R. mucilaginosa tolerance to such a broad range of Pb²⁺ concentrations can be attributed to energy adaption. In contrast, high Pb²⁺ stress caused ATP deficiency. Then, the subsequent degradation of intracellular defense systems further intensified Pb toxicity. This study correlated responses of EMP, TCA, and OXPHOS pathways in R. mucilaginosa under Pb stress, hence providing new insights into the fungal resistance to heavy metal stress. IMPORTANCE Glycolysis (EMP), the tricarboxylic acid (TCA) cycle, and oxidative phosphorylation (OXPHOS) are critical metabolism pathways for microorganisms to obtain energy during the resistance to heavy metal (HM) stress. However, these pathways at the genetic level have not been elucidated to evaluate their cytoprotective functions for Rhodotorula mucilaginosa under Pb stress. In this study, we investigated these three pathways based on biochemical analysis and RNA sequencing. Under low-medium (500-1,000 mg · L^-1) Pb²⁺ stress, ATP production was stimulated mainly due to the upregulation of genes associated with the TCA cycle and the electron transport chain (ETC). Such an energy compensatory mechanism could allow R. mucilaginosa acclimation to a broad range of Pb²⁺ concentrations (up to 1000 mg · L^-1). In contrast, high (2500 mg · L^-1) Pb²⁺ stress exerted its excessive toxicity by provoking ATP deficiency and damage to intracellular resistance systems. This study provided new insights into R. mucilaginosa resistance to HM stress from the perspective of metabolism.

Single-cell RNA sequencing of freshly isolated bovine milk cells and cultured primary mammary epithelial cells

Bovine mammary function at molecular level is often studied using mammary tissue or primary bovine mammary epithelial cells (pbMECs). However, bulk tissue and primary cells are heterogeneous with respect to cell populations, adding further transcriptional variation in addition to genetic background. Thus, understanding of the variation in gene expression profiles of cell populations and their effect on function are limited. To investigate the mononuclear cell composition in bovine milk, we analyzed a single-cell suspension from a milk sample. Additionally, we harvested cultured pbMECs to characterize gene expression in a homogeneous cell population. Using the Drop-seq technology, we generated single-cell RNA datasets of somatic milk cells and pbMECs. The final datasets after quality control filtering contained 7,119 and 10,549 cells, respectively. The pbMECs formed 14 indefinite clusters displaying intrapopulation heterogeneity, whereas the milk cells formed 14 more distinct clusters. Our datasets constitute a molecular cell atlas that provides a basis for future studies of milk cell composition and gene expression, and could serve as reference datasets for milk cell analysis.

Extending the lymphoblastoid cell line model for drug combination pharmacogenomics

Combination drug therapies have become an integral part of precision oncology, and while evidence of clinical effectiveness continues to grow, the underlying mechanisms supporting synergy are poorly understood. Immortalized human lymphoblastoid cell lines (LCLs) have been proven as a particularly useful, scalable and low-cost model in pharmacogenetics research, and are suitable for elucidating the molecular mechanisms of synergistic combination therapies. In this review, we cover the advantages of LCLs in synergy pharmacogenomics and consider recent studies providing initial evidence of the utility of LCLs in synergy research. We also discuss several opportunities for LCL-based systems to address gaps in the research through the expansion of testing regimens, assessment of new drug classes and higher-order combinations, and utilization of integrated omics technologies.

The Transcriptome of SH-SY5Y at Single-Cell Resolution: A CITE-Seq Data Analysis Workflow

Cellular Indexing of Transcriptomes and Epitopes by Sequencing (CITE-seq) is a recently established multimodal single cell analysis technique combining the immunophenotyping capabilities of antibody labeling and cell sorting with the resolution of single-cell RNA sequencing (scRNA-seq). By simply adding a 12-bp nucleotide barcode to antibodies (cell hashing), CITE-seq can be used to sequence antibody-bound tags alongside the cellular mRNA, thus reducing costs of scRNA-seq by performing it at the same time on multiple barcoded samples in a single run. Here, we illustrate an ideal CITE-seq data analysis workflow by characterizing the transcriptome of SH-SY5Y neuroblastoma cell line, a widely used model to study neuronal function and differentiation. We obtained transcriptomes from a total of 2879 single cells, measuring an average of 1600 genes/cell. Along with standard scRNA-seq data handling procedures, such as quality checks and cell filtering procedures, we performed exploratory analyses to identify most stable genes to be possibly used as reference housekeeping genes in qPCR experiments. We also illustrate how to use some popular R packages to investigate cell heterogeneity in scRNA-seq data, namely Seurat, Monocle, and slalom. Both the CITE-seq dataset and the code used to analyze it are freely shared and fully reusable for future research.

Single-cell RNA-seq reveals transcriptomic heterogeneity mediated by host-pathogen dynamics in lymphoblastoid cell lines

Lymphoblastoid cell lines (LCLs) are generated by transforming primary B cells with Epstein–Barr virus (EBV) and are used extensively as model systems in viral oncology, immunology, and human genetics research. In this study, we characterized single-cell transcriptomic profiles of five LCLs and present a simple discrete-time simulation to explore the influence of stochasticity on LCL clonal evolution. Single-cell RNA sequencing (scRNA-seq) revealed substantial phenotypic heterogeneity within and across LCLs with respect to immunoglobulin isotype; virus-modulated host pathways involved in survival, activation, and differentiation; viral replication state; and oxidative stress. This heterogeneity is likely attributable to intrinsic variance in primary B cells and host–pathogen dynamics. Stochastic simulations demonstrate that initial primary cell heterogeneity, random sampling, time in culture, and even mild differences in phenotype-specific fitness can contribute substantially to dynamic diversity in populations of nominally clonal cells.

Single-Cell Expression Variability Implies Cell Function

As single-cell RNA sequencing (scRNA-seq) data becomes widely available, cell-to-cell variability in gene expression, or single-cell expression variability (scEV), has been increasingly appreciated. However, it remains unclear whether this variability is functionally important and, if so, what are its implications for multi-cellular organisms. Here, we analyzed multiple scRNA-seq data sets from lymphoblastoid cell lines (LCLs), lung airway epithelial cells (LAECs), and dermal fibroblasts (DFs) and, for each cell type, selected a group of homogenous cells with highly similar expression profiles. We estimated the scEV levels for genes after correcting the mean-variance dependency in that data and identified 465, 466, and 364 highly variable genes (HVGs) in LCLs, LAECs, and DFs, respectively. Functions of these HVGs were found to be enriched with those biological processes precisely relevant to the corresponding cell type’s function, from which the scRNA-seq data used to identify HVGs were generated—e.g., cytokine signaling pathways were enriched in HVGs identified in LCLs, collagen formation in LAECs, and keratinization in DFs. We repeated the same analysis with scRNA-seq data from induced pluripotent stem cells (iPSCs) and identified only 79 HVGs with no statistically significant enriched functions; the overall scEV in iPSCs was of negligible magnitude. Our results support the “variation is function” hypothesis, arguing that scEV is required for cell type-specific, higher-level system function. Thus, quantifying and characterizing scEV are of importance for our understating of normal and pathological cellular processes.