Graphical Abstract

   graphic file with name fx1.jpg
   [35]Open in a new tab

Highlights

     * •
       scRNA-seq reveals that nicotine alters transcriptome of hESC
       differentiation
     * •
       Nicotine reduces cell viability, increases ROS level, and changes
       cell cycle of EBs
     * •
       Nicotine disturbs intracellular Ca^2+ in hESC-derived
       cardiomyocytes
     __________________________________________________________________

   Using single-cell RNA sequencing (scRNA-seq), Wu and colleagues
   elucidate the adverse effects of nicotine exposure on various cell
   lineages derived from human embryonic stem cells. This study offers new
   insights into the toxicity of environmental factors on human embryonic
   development.

Introduction

   Maternal smoking during pregnancy is an established risk factor for
   birth defects such as miscarriage, growth restriction, and premature
   birth ([36]Jaddoe et al., 2008). It is closely associated with adverse
   neurobehavioral, cardiovascular, respiratory, endocrine, and metabolic
   outcomes in the offspring, which can persist into adulthood
   ([37]Holbrook, 2016). Nicotine, the main chemical constituent of
   tobacco smoking, is primarily responsible for the elevated risk
   ([38]Holbrook, 2016). Unfortunately, the introduction and spread of new
   tobacco products containing nicotine, such as e-cigarettes, is
   reversing recent progress toward reduction of tobacco use ([39]Bao
   et al., 2018).

   A large body of research has elucidated the negative effects of
   nicotine in animals, mainly in rodent models. Animal studies have
   demonstrated that nicotine exposure during pregnancy has detrimental
   effects on fetal development, such as cellular damage, increased
   inflammation ([40]Mohsenzadeh et al., 2014), oxidative stress ([41]Lin
   et al., 2014), endoplasmic reticulum stress ([42]Wong et al., 2016),
   and impaired cell replication ([43]Repo et al., 2014, [44]Slotkin
   et al., 1987). The suitability of clinical translation of these
   studies, however, remains questionable due to interspecies
   physiological differences and uncertainty over the degree and route of
   nicotine exposure ([45]Tizabi, 2007, [46]Winzer-Serhan, 2008). To
   address these issues, some studies have attempted to study the effects
   of nicotine using human cells. For example, using microarray analysis,
   [47]Liszewski et al. (2012) demonstrated that tobacco smoke and
   nicotine have lineage- and stage-specific effects on differentiated
   human embryonic stem cell (hESCs).

   Although the in vitro differentiation of embryonic body (EB) model can
   be used to mimic early developments from pre-implantation epiblasts to
   lineage-committed progenitors, conventional bulk RNA sequencing
   (RNA-seq) analysis has limitations for studying the individual cellular
   heterogeneity within the EBs. With the recent advent of
   microdroplet-based single-cell RNA-seq (scRNA-seq) technologies, it is
   now feasible to analyze transcriptomes at the single-cell level within
   heterogeneous cell populations ([48]Blakeley et al., 2017, [49]Paik
   et al., 2018). Here, we used scRNA-seq of EBs to characterize the
   effects of nicotine on hESC differentiation. We found that nicotine
   exposure reduced cell viability and increased reactive oxygen species
   (ROS), resulting in aberrant formation and differentiation of EBs.
   Nicotine exposure also altered cell cycling in endothelial, stromal,
   and muscle progenitor cells differentiated from hESCs. Furthermore,
   nicotine caused lineage-specific effects and dysregulated cell-to-cell
   communication. We found abnormal Ca^2+ signaling pathways in muscle
   cells upon nicotine exposure that was verified using hESC-derived
   cardiomyocytes. Taken together, the effects of nicotine exposure on
   hESC differentiation at the single-cell transcriptomic level offer new
   insights into mechanisms of nicotine toxicity on early embryonic
   development, and can provide new tools for optimizing drug toxicity
   screening.

Results

scRNA-Seq Analysis Reveals Six Major Types of Progenitor Cells

   To investigate the effects of nicotine on hESC differentiation, we
   performed microdroplet-based scRNA-seq to identify unique cell lineages
   on day 21 control and nicotine-exposed EBs ([50]Figure 1A). We used
   10 μM nicotine exposure for 21 days, which is similar to nicotine
   concentrations found in fetal serum ([51]Luck et al., 1985) and has
   been used in prior hESC studies ([52]Hirata et al., 2016,
   [53]Zdravkovic et al., 2008). After dissociation, transcriptomic data
   of 5,646 single cells from nicotine-exposed EBs and 6,847 single cells
   from control EBs were acquired. Sequenced data showed high read depth,
   and were mapped to approximately 3,000 median genes per cell
   ([54]Figure S1A, left). The percentage of mitochondrial genes present
   in most cells was less than 10% ([55]Figure S1A, right). We used the
   Seurat package ([56]Satija et al., 2015) to perform principal-component
   analysis and t-distributed stochastic neighbor embedding (t-SNE)
   analysis. Control EBs were divided into 13 clusters, and
   nicotine-exposed EBs were divided into 12 clusters that exhibited
   distinct gene expression patterns ([57]Figures S1B and S1C). Control
   and nicotine-exposed EBs contained similar cell-type markers, without
   any observed differences in cell types between the two samples
   ([58]Figure S1B).

Figure 1.

   [59]Figure 1
   [60]Open in a new tab

   scRNA-Seq Analysis Reveals Cell Lineages in Control and
   Nicotine-Exposed Embryoid Bodies

   (A) Process flow diagram of scRNA-seq analysis on hESC differentiation.
   Single cells were collected from two independent EB differentiation
   experiments from day 21 EBs (nicotine-exposed versus control) and were
   prepared by single-cell barcoded droplets and chemicals from 10×
   Genomics. Bioinformatics data were processed using Seurat. Cell-type
   marker, differentially expressed gene, cell communication, and pathway
   analyses were performed to investigate the effects of nicotine exposure
   on hESC differentiation.

   (B) Separated (left) and combined (middle and right) t-SNE plots of
   single cells from control and nicotine-exposed EBs. We defined six main
   types of progenitor cells in day 21 EBs, including muscle progenitor
   cells (clusters 3 and 13), liver progenitor cells (cluster 5), neural
   progenitor cells (clusters 3, 4, 8, and 11), stromal progenitor cells
   (cluster 6), epithelial progenitor cells (clusters 2 and 12), and
   endothelial progenitor cells (cluster 11). In addition,
   undifferentiated stem-like cells (USCs) (cluster 1) and undetermined
   cells (UDCs) (cluster 9) were also identified.

   (C) Heatmap showing the expression pattern of top 10 differential genes
   in each cell type. Representative differential genes for each cell type
   are listed on the right side. The complete lists of differential genes
   for each cell type are listed in [61]Table S3.

   (D) Violin plots show the expression level distributions of marker
   genes across cell types.

   Next, we performed integrative analysis to compare the cell proportions
   and gene expression differences in each cell type between nicotine and
   control EBs. Nicotine exposure induced widespread transcriptomic
   changes, which were manifested as a shift in the t-SNE projections of
   singlets ([62]Figure 1B, middle). Previous reports with bulk RNA-seq
   data also indicated that nicotine affects gene expression in multiple
   cell lineages ([63]Liszewski et al., 2012). A total of 13 individual
   clusters were defined from the combined datasets (termed C1 to C13)
   ([64]Figure 1B, right). Based on differential genes enriched in each
   cluster, six major types of progenitor cells were identified by Seurat
   ([65]Figures 1C and 1D). Clusters 3, 4, 8, and 10 were associated with
   high expression of LHX2 and NR2F1 ([66]de Melo et al., 2016, [67]Tang
   et al., 2010) and annotated as neural cells. Cluster 5 represented
   liver progenitor cells with a high expression of FRZB and PTN
   ([68]Michelotti et al., 2014, [69]Shen et al., 2015). Cluster 6 was
   annotated as stromal progenitor cells with a high expression of SFRP2
   and COL2A1 ([70]Saito et al., 2013, [71]Tabib et al., 2018). Cluster 11
   was annotated as endothelial progenitor cells with high expression of
   GDF15 and DDIT3 ([72]Ahrens et al., 2011, [73]Loinard et al., 2012).
   Cluster 2 and 12 were annotated as epithelial progenitor cells (EpiPCs)
   with high expression of IGFBP5 and PODXL ([74]Sugrue et al., 2016,
   [75]Zhu et al., 2016). Clusters 7 and 13 showed a high expression of
   HAPLN1 and S100A11, and were annotated as muscle progenitor cells
   ([76]DeLaughter et al., 2013, [77]Malmstrom et al., 2004). Cluster 1
   was enriched for pluripotency genes such as TERF1 and POU5F1, and was
   annotated as “undifferentiated stem-like cells”. Cluster 9 was enriched
   for cytoskeletal genes such as ACTB and TUBB, and was annotated as
   “undetermined cells”.

   To further confirm our cluster annotation, we found genes specifically
   expressed in each cell type that were enriched for the expected
   appropriate gene ontology (GO) terms. For example, genes that were
   specifically expressed in muscle progenitor cell clusters were
   significantly enriched for the cytosolic Ca^2+ pathway (p = 2.92 ×
   10^−11) and skeletal system development (p = 3.36 × 10^−4).
   Genes expression in the neural progenitor cell clusters were
   significantly enriched for nervous system development (p = 1.01 ×
   10^−6) and sensory organ development (p = 5.60 × 10^−8). Genes
   expression in the liver progenitor cell cluster were enriched for liver
   development (p = 2.52 × 10^−4) and response to lipid (p = 1.74 ×
   10^−2). Genes expression in the endothelial progenitor cell cluster
   were significantly enriched for blood vessel development (p = 1.3 ×
   10^−4) and angiogenesis (p = 1.61 × 10^−3). Genes expressions in the
   EpiPC clusters were enriched for lung development (p = 2.43 × 10^−5)
   and kidney development (p = 5.67 × 10^−4) ([78]Figure S1D).

   It should be noted that neural, muscle, and epithelial progenitor cells
   consisted of several sub-clusters. Neural progenitor cells were further
   divided into four subsets: clusters 3, 4, 8, and 10. Cluster 3 showed a
   high expression of LHX5/HESX1 that is related to forebrain development
   ([79]Martynova et al., 2018, [80]Zhao et al., 1999). Cluster 4 was
   enriched for HMGB2 and PTTG1, which are highly expressed in
   proliferating neural stem cells ([81]Kimura et al., 2018). Cluster 8
   was annotated as neural progenitor cells with an enrichment of HNRNPH1
   and PTPRS ([82]Tchetchelnitski et al., 2014, [83]Yazdani et al., 2015),
   which are related to sensory neurons development. Cluster 10 was
   enriched for LHX2 and NR2F1 and expressed eye development genes ([84]de
   Melo et al., 2016, [85]Tang et al., 2010). Muscle progenitor cells were
   divided into two subsets: clusters 7 and 13. Cluster 7 was annotated as
   muscle cells for the expression of HAPLN1 and S100A11, which are highly
   expressed in smooth muscle cells ([86]DeLaughter et al., 2013,
   [87]Malmstrom et al., 2004). Cluster 13 was enriched for ZFHX3 and
   NR2F2, which are related to cardiac muscle development ([88]Berry
   et al., 2001, [89]Pei et al., 2017). EpiPCs were divided into two
   subsets: clusters 2 and 12. Cluster 2 was enriched for IGFBP5 and HES1,
   which are related to eye development ([90]Liu et al., 2013, [91]Sugrue
   et al., 2016). Cluster 12 was enriched for B3GNT7 and PODXL, which are
   highly expressed in stem-like epithelial cells ([92]Dumont-Lagace
   et al., 2017) ([93]Figures 1D and [94]S1D).

   Overall, six major types of progenitors (neural, liver, stromal,
   endothelial, epithelial, and muscle) were identified from scRNA-seq
   data of EBs based on cell markers detected by Seurat. These data may be
   useful for modeling nicotine exposure on individual organs and cells
   within the developing fetus.

Nicotine Elicits Cell-Type-Specific Response in Differentiated EBs

   Integrated analysis of control and nicotine-exposed EBs at the
   single-cell level enables us to quantitatively assess
   cell-type-specific responses to nicotine. Quantification of the
   cell-type compositional changes showed changes from 5% reduction in
   epithelial progentior cells to 4% increase in liver progenitor cells
   following nicotine exposure ([95]Figures 2A and [96]S2A). Next, we
   performed comparative analysis and calculated the average expression of
   both the nicotine-exposed and control cells to determine differentially
   expressed genes (DEGs) in each cell type ([97]Figure 2B).
   Interestingly, there was a marked difference in the number of DEGs
   among different cell types, ranging from 5 to 103 genes with a p-value
   less than 0.01 and a log fold-change more than 0.25. For example, we
   observed 5 DEGs in liver progenitor cells with nicotine treatment,
   whereas muscle progenitor cells exhibited the greatest number of 103
   DEGs. Among these DEGs, BNIP3, and metallothionein family genes (MT1X,
   MT1G, MT1E, and MT2A) were uniformly downregulated and RPS10 was
   upregualted in most of cell types ([98]Figure 2B). BNIP3 gene is an
   important regulator during long-term nicotine-induced cell death in
   several cell types ([99]Erkan et al., 2005, [100]Tang et al., 2007).
   Metallothionein family genes play a role in the protection against
   metal toxicity and oxidative stress, and have been shown to be
   suppressed in chronic smokers ([101]Billatos et al., 2018). These genes
   are involved in apoptosis, ROS generation, mitochondrial function, and
   response to metal ion pathways ([102]Figure 2C), indicating that EBs
   have poor cell survival upon nicotine exposure.

Figure 2.

   [103]Figure 2
   [104]Open in a new tab

   Nicotine Exposure Induces Cell-Type-Specific Response

   (A) Cell proportion fluctuation for each cell type with nicotine
   exposure. Cell proportion fluctuation for each cluster with nicotine
   exposure is shown in [105]Figure S2A. USCs, undifferentiated stem-like
   cells; UDCs, undetermined cells; EpiPCs, epithelial progenitor cells;
   ECs, endothelial progenitor cells.

   (B) Plots of the average expression of genes from control and
   nicotine-exposed EBs for each cell type. Significant differentially
   expressed genes are labeled in the plots (p < 0.05). The complete list
   of differentially expressed genes upon nicotine exposure for each
   cluster is shown in [106]Table S5.

   (C) Pathway enrichment analysis of statistically significant gene
   ontologies following nicotine exposure. The size of the circle
   represents the significance of gene ontologies, and the darkness of
   color represents the number of genes involved in the gene ontologies.
   The complete lists of differentially expressed gene-related pathways
   upon nicotine exposure for each cluster are listed in [107]Table S6.

   (D) Heatmap for differentially expressed genes in scRNA-seq data,
   public hepatic cell line (HepaRG), human gingival epithelium cell line
   (HGEC), human smooth muscle cells (HSMCs), and human iPSC-derived
   endothelial cell (hiPSC-ECs) in terms of fold-change (nicotine-exposed
   relative to control EBs). Each row represents a single differentially
   expressed gene identified in (B). The names of cell types are labeled
   on the left side. Differentially expressed genes in scRNA-seq data
   corresponding to the public datasets are labeled with a rectangle.

   See also [108]Figure S2.

   Nicotine also showed cell-type-specific responses. APOE, TUBA1A, and
   NDUFC2 were significantly upregulated, and H1F0 and SRRM2 were
   downregulated, in neural progenitor cells ([109]Figures 2B and 2C).
   Abnormal expression of these genes can lead to β-amyloid formation and
   increased synaptic transmission ([110]Moreno-Gonzalez et al., 2013),
   brain malformations ([111]Aiken et al., 2017), and intellectual
   disability ([112]Tanaka et al., 2018). In muscle progenitor cells, the
   most upregulated gene following nicotine exposure was HSP90AA1, a
   myosin chaperone protein gene involved in muscle development and
   disease ([113]Armant et al., 2016, [114]Etard et al., 2015). Increased
   expression of HMGB1, known to regulate cardiac excitation-contraction
   coupling by enhancing the sarcoplasmic reticulum Ca^2+ leakage through
   Toll-like receptor 4 (TLR4)-ROS signaling in cardiac muscle cells, was
   also observed ([115]Figures 2B, 2C, and [116]S2C).

   In stromal progenitor cells, LDHA and DAPL1, known to regulate nutrient
   levels and amino acid acetylation, were downregulated upon nicotine
   exposure ([117]Figures 2B and 2C). BINP3, related to lipid metabolism,
   was downregulated in liver cells ([118]Glick et al., 2012). In
   epithelial progenitor cells, macrophage migration inhibitory factor
   (MIF), which is associated with chronic obstructive pulmonary disease
   in human, was downregulated ([119]Sauler et al., 2015). WDR77, required
   for proliferation of lung and prostate epithelial cells during
   development and tumorigenesis ([120]Sheng and Wang, 2016), was
   upregulated in epithelial progenitor cells in the presence of nicotine
   ([121]Figures 2B and 2C). In endothelial progenitor cells, LDHA, DDIT3,
   and IFITM1 were downregulated, and HNRNP2B1 was upregulated upon
   nicotine exposure ([122]Figures 2B and 2C). Downregulation of LDHA is
   related to the suppression of glycolysis and endothelial cell
   dysfunction ([123]Xu et al., 2016). Downregulated DDIT3 exhibits
   reduced ER stress response upon long-term cigarette smoke exposure
   ([124]Geraghty et al., 2011). IFITM1 downregulation is related to
   endothelial lumen formation during angiogenesis ([125]Popson et al.,
   2014) ([126]Figures 2B and 2C). Likewise, we confirmed
   cell-type-specific responses to long-term nicotine exposure on EBs
   using GO pathway enrichment analysis ([127]Figure 2C).

   To determine how our analysis correlates with previously reported bulk
   RNA-seq data from specific cell lines, we examined the expression
   fold-changes of DEGs in published gene expression data after nicotine
   exposure, including hepatic cells (HepaRG), human gingival epithelium
   cells (HGECs), human smooth muscle cells (HSMCs), and human
   iPSC-derived endothelial cells (hiPSC-ECs) ([128]De Abrew et al., 2016,
   [129]Gumus et al., 2008, [130]Yoshiyama et al., 2014). Interestingly,
   the trend in gene expression changes in cluster12 epithelial progenitor
   cell was similar to that of HGECs, whereas an opposite trend of gene
   expression changes was observed between cluster-2 epithelial progenitor
   cells and HGECs. This may be due to HGECs being more similar to
   stem-like “cluster-12 epithelial cell” in biological identity and more
   different from “cluster-2 epithelial cell”. In HSMCs, the overall
   change of DEGs was subtle, but the trends in fold-change expression
   were consistent with EB-derived muscle cells. The genes downregulated
   in EB-derived endothelial cells also reduced their expression in
   hiPSC-ECs upon nicotine exposure ([131]Figure 2D). Taken together, DEG
   analysis showed cell-type-specific transcriptomic changes upon nicotine
   exposure, which are consistent with previously reported bulk RNA-seq or
   microarray analysis in different cell types ([132]De Abrew et al.,
   2016, [133]Gumus et al., 2008, [134]Yoshiyama et al., 2014). Our data
   thus provide a novel method for evaluating nicotine toxicity in
   heterogeneous populations of human EBs at a single-cell level.

Nicotine Dysregulates Viability, ROS Generation, and Cell Cycle in EBs

   Infants exposed to nicotine prenatally often exhibit lower birth
   weights than their peers ([135]Fried and Oconnell, 1987, [136]Slotkin,
   1998). Animal studies have shown that nicotine exposure during
   pregnancy induces cellular damage, oxidative stress, and impaired cell
   replication ([137]Repo et al., 2014, [138]Slotkin et al., 1987).
   However, the molecular mechanisms remain poorly understood. Our DEG
   analysis showed that long-term nicotine exposure induced apoptosis and
   ROS generation mediated by the downregulation of BNIP3 and
   metallothionein family genes ([139]Figures 2B and 2C). Therefore, we
   performed several assays to confirm decreased survival of
   nicotine-exposed EB. Nicotine-exposed EBs were smaller than control EBs
   ([140]Figures 3A and 3B), and cell viability was significantly reduced
   based on quantification of ATP, an indicator of metabolically active
   cells ([141]Figure 3C). We also found higher levels of ROS in
   nicotine-exposed EBs compared with control EBs ([142]Figure 3D).

Figure 3.

   [143]Figure 3
   [144]Open in a new tab

   Nicotine Reduces Cell Viability, Increases ROS Levels, and Changes Cell
   Cycle in EBs

   (A) Representative bright-field images of day 21 control and
   nicotine-exposed EBs. Scale bar, 100 μm.

   (B) Size measurement of control (n = 45) and nicotine-exposed (n = 43)
   EBs in terms of diameter. EBs were collected from three independent EB
   differentiation experiments and pooled together for size measurement.
   ^∗∗∗p < 0.001.

   (C) Cell viability assay of EBs based on quantitation of the ATP
   present in control and nicotine-exposed EBs. Cell viability were
   measured from three independent experiments. ^∗p < 0.05.

   (D) Reactive oxygen species (ROS) generation in control and
   nicotine-exposed EBs. ROS was measured from three independent
   experiments. ^∗p < 0.05.

   (E) Proportion of cluster 13 muscle cells in G2M, S, or G1 phase.

   (F) Distribution of TUBB4B expression in three cell-cycle phases in
   control (top) and the nicotine-exposed EBs (bottom). Dashed lines
   represent the center of G2M phase in control and nicotine-exposed EBs.

   Clinical and animal studies have shown that nicotine exposure changes
   the dynamics of cell replication and causes growth restriction
   ([145]Repo et al., 2014). We thus analyzed cell cycling in the
   scRNA-seq data to evaluate the growth of EB after nicotine treatment by
   calculating cell-cycle phase scores based on canonical markers
   ([146]Nestorowa et al., 2016). Relative to control EBs,
   nicotine-exposed EBs exhibited a 12% decrease in G1 phase, a 6%
   increase in G2M phase, and a 5.5% increase of S phase in endothelial
   progenitor cells. In stromal progenitor cells, we found an 11% decrease
   in G1 phase and a 12% increase in S phase ([147]Figure S2B).
   Surprisingly, we found that there was a 20% decrease of cells in the G1
   phase, a 5% increase of cells in the G2M phase, and a 15% increase of
   cells in S phase in the “cluster-13 muscle progenitor cell”
   ([148]Figure 3E). For example, TUBB4B, a G2M phase marker, was
   differentially expressed in muscle progenitor cells from
   nicotine-exposed EBs versus control EBs ([149]Figure 3F). Consequently,
   nicotine exposure increased ROS production and cell death in EBs and
   affected the cell cycle of endothelial, stromal, and muscle progenitor
   cells.

Nicotine Exposure Dysregulates Cell-to-Cell Communication of Differentiated
EBs

   Smoking and nicotine consumption increase the pathological risk in
   endocrine, reproductive, respiratory, cardiovascular, and neurologic
   systems that all rely on intricate and dynamic interactions among
   multiple functional cell types for homeostasis and function
   ([150]Kawasaki et al., 2011, [151]Rehan et al., 2009). The effect of
   nicotine on cell-to-cell communication, however, is not well
   understood. Recent studies using in vitro co-cultured systems indicate
   that cell-to-cell communication could be affected by nicotine exposure
   ([152]Holownia et al., 2015, [153]Larsen et al., 2016, [154]Liu et al.,
   2017). Our study used a dataset of human ligand-receptor pairs
   ([155]Ramilowski et al., 2015) to define intercellular communication
   networks. To examine the effects of nicotine on cell-to-cell
   communication, we also analyzed ligand-receptor expression differences
   in nicotine-exposed EBs and control EBs. Overall, we observed increased
   intercellular communication for each EB cell type upon nicotine
   exposure ([156]Figures 4A and 4B). For example, the number of
   ligand-receptor pairs in autocrine circuits from muscle progenitor cell
   was increased from 51 to 85, and the number of ligand-receptor pairs in
   muscle-neuron crosstalk was increased from 47 to 71 ([157]Figure 4B).

Figure 4.

   [158]Figure 4
   [159]Open in a new tab

   Nicotine Increased Cell-to-Cell Communication and Induced HMGB1-TLR4
   Pathway among Cell Types in EBs

   (A) Intercellular communication analysis among cell types in the
   control and nicotine-exposed EBs. Line color indicates ligands
   broadcast by the cell population of the same color (labeled). Lines
   connect to cell populations where cognate receptors are expressed. Line
   thickness is proportional to the number of ligands where cognate
   receptors are present in the recipient cell population. Loops indicate
   autocrine circuits. Map quantifies potential communication but does not
   account for anatomic position or boundaries of cell populations.
   Crosstalk among cell types in control and nicotine-exposed EBs are
   listed in [160]Table S7.

   (B) Detailed view of ligands broadcast by each cell type and those
   populations expressing cognate receptors primed to receive a signal.
   Numbers indicate the quantity of ligand-receptor pairs for each
   inter-population link. Colors of the cell population corresponds to (A)
   above.

   (C) Venn diagrams of ligands present in control or nicotine-exposed EBs
   and differentially expressed genes between control and nicotine-exposed
   EBs.

   (D) Heatmap of the expression of differentially expressed ligands in
   control versus nicotine-exposed EBs across different cell types in
   terms of average fold-change of expression. Seven ligands present in
   both control and nicotine-exposed EBs (labeled “both”) were selected as
   differentially expressed ligands. Three ligands were selected as
   expressed only in nicotine-exposed EBs and labeled as “nicotine.”

   (E) Intercellular communication analysis among different cell types
   shows activated HMGB1-TLR4 pathway upon nicotine exposure. The meaning
   of the thickness and color of line is explained in (A). USCs,
   undifferentiated stem-like cells; UDCs, undetermined cells; EpiPCs,
   epithelial progenitor cells; ECs, endothelial progenitor cells.

   Next, we analyzed the expression of ligands in each cell type to
   identify 61 ligands that were expressed in both nicotine-exposed EBs
   and control EBs, of which 7 were differentially expressed. Seventy-five
   ligands were mainly expressed in nicotine-exposed EBs, of which 3 were
   enriched in certain cell types identified by DEG analysis ([161]Figures
   4C and 4D). One extracellular ligand, high-mobility group box 1
   (HMGB1), was uniformly upregulated in multiple cell types of
   nicotine-exposed EBs, but was not presented in control EBs
   ([162]Figures 4C and 4D). Previous research has shown that HMGB1
   regulates Ca^2+ handling and cellular contractility by activating its
   receptor Toll-like receptor 4 (TLR4) in rat cardiomyocytes, which plays
   an important role in the pathogenesis of cardiac dysfunction in many
   diseases ([163]Zhang et al., 2014). Here, we found that HMGB1-TLR4
   signaling, although not specific, was also activated in muscle
   progenitor cells ([164]Figure 4E), suggesting that the activated
   HMGB1-TLR4 pathway may play an important role in cardiac dysfunction
   upon nicotine exposure.

Disturbance of Intracellular Ca^2+ Handling in hESC-Derived Cardiomyocytes by
Nicotine Exposure

   Our data indicate that nicotine affects the expression of genes
   associated with intracellular Ca^2+ handling via the HMBG1-TLR4 pathway
   in cardiac muscle cells, and animal studies have shown that nicotine
   exposure disrupts intracellular Ca^2+ homeostasis in cardiac cells
   ([165]Hu et al., 2013). To investigate whether nicotine affects the
   Ca^2+ handling in cardiac muscle cells, we first checked the expression
   of HMGB1 and TLR4 in hESC-derived cardiomyocytes. The expression of
   HMGB1 was increased by 2-fold, and TLR4 increased by 90-fold in
   hESC-derived cardiomyocytes exposed to nicotine ([166]Figure 5A). We
   next conducted single-cell Ca^2+ measurement using Fura-2 in
   hESC-derived cardiomyocytes ([167]Figure S3). As shown in [168]Figures
   5B–5E, nicotine increased the diastolic Ca^2+ ([169]Figures 5B and 5C)
   and reduced the Ca^2+ transient amplitude ([170]Figure 5D), accompanied
   by prolonged Ca^2+ decay ([171]Figure 5E), suggesting compromised
   intracellular Ca^2+ homeostasis. We found that nicotine increased the
   propensity for arrhythmic Ca^2+ release in hESC-derived cardiomyocytes,
   as indicated by the arrows in [172]Figures 5F and 5G. These data
   strongly suggest that nicotine increases Ca^2+-release abnormalities at
   the cellular level, predisposing these cells to Ca^2+-associated
   arrhythmia.

Figure 5.

   [173]Figure 5
   [174]Open in a new tab

   Nicotine Disrupts Intracellular Ca^2+ Handling in hESC-Derived
   Cardiomyocytes

   (A) Real-time PCR verification of the expression levels of HMGB1 and
   TLR4 in hESC-derived cardiomyocytes for nicotine-exposed and control
   group. Real-time PCR was conducted from three independent
   differentiation experiments. Nicotine and ethanol were added during
   cardiomyocyte differentiation process. Ethanol-exposed hESC-derived
   cardiomyocytes were used as a control. ^∗∗∗p < 0.001.

   (B–E) Representative single-cell Ca^2+ traces with and without
   treatment of nicotine (1 μM) (B). Cells were paced at 0.5 Hz and
   recorded at 37°C. Nicotine affects diastolic Ca^2+ (C), Ca^2+ transient
   amplitude (D), and Ca^2+ uptake (E).

   (F) Representative arrhythmic Ca^2+ release. The arrhythmic Ca^2+
   release is indicated by the arrows.

   (G) Percentage of cells with normal or arrhythmic Ca^2+ release.
   Numbers indicate the recorded cell number in each condition. Recording
   was conducted with two independent differentiation. The cell number
   shows the pooled data.

   ^∗∗p < 0.01, ^∗∗∗p < 0.001. Error bars indicate SEM.

Discussion

   In this study, we performed scRNA-seq analysis on a total of 12,500
   single cells generated from human ESC-derived EBs following 21 days of
   culture with or without nicotine ([175]Figure 1). Previous studies have
   demonstrated that nicotine concentrations in fetal serum are much
   higher than in maternal serum, ranging from 0.3 to 15.4 μM in fetal
   serum ([176]Luck et al., 1985). Based on these reports, and the
   concentrations studied in other investigations of nicotine effect on
   hESCs (0.1–10.0 μM), we decided to use 10 μM of nicotine during hESC
   differentiation ([177]Hirata et al., 2016, [178]Zdravkovic et al.,
   2008). We also found that a 6-day exposure to nicotine reduces
   viability in hESCs ([179]Figure S3D), suggesting that nicotine affects
   embryo development as early as the pre-implantation stage.

   We did not observe cell-type differences between nicotine-exposed EBs
   and control EBs, although there were minor changes in the cell-type
   distribution upon nicotine exposure. However, DEG patterns from various
   progenitor cell populations indicated broad effects on cells derived
   from all three germ layers (neural, stromal, muscle, endothelial, and
   epithelial progenitor cells). This is consistent with clinical
   observations that nicotine-exposed infants have health problems
   throughout their lives, including impaired function of the endocrine,
   reproductive, respiratory, cardiovascular, and neurologic systems
   ([180]Warren et al., 2014). In addition, although the current
   technology does not allow us to conduct proteomic analyses at a
   single-cell resolution level, we searched available nicotine-associated
   proteomic datasets in PubMed and found that our DEGs, HSPA8, BAX, and
   CKB were also reported to be upregulated in mouse neurons
   ([181]Matsuura et al., 2016), MIF was downregulated in human
   endothelial cells ([182]Zhang et al., 2014), and that SLC25A4, REEP5,
   and ATP5F1 were upregulated in human epithelial cells ([183]Ghosh
   et al., 2018) ([184]Figure S2E). Moreover, among these DEGs, BNIP3 is
   uniformly downregulated in most cell types ([185]Figure 2B). We
   observed downregulated BNIP3 expression in nicotine-exposed EBs
   compared with the control EBs at the protein level ([186]Figure S2D).
   These lines of evidence suggest that our scRNA-seq data can reflect
   protein expression to some extent.

   Our scRNA-seq analysis suggests that downregulation of BNIP3 and
   metallothionein family gene expression in multiple progenitor cell
   types upon nicotine exposure provides the molecular mechanisms that
   lead to altered DEG patterns. It has been reported that reduced BNIP3
   expression correlates with poor cell survival following long-term
   nicotine exposure ([187]Tang et al., 2007). Our data indicate that
   hESC-derived cells adapt to long-term nicotine exposure and cell damage
   resulting from downregulation of BNIP3 ([188]Figure 2B). This may be
   analogous to the observed low birth weight, preterm birth, and
   perinatal death in a developing fetus as a result of maternal smoking
   during pregnancy ([189]Warren et al., 2014).

   Nicotine also displayed cell-type-specific adverse effects, consistent
   with previous findings in animal and clinical studies ([190]Holbrook,
   2016). For example, we found that APOE was upregulated in
   nicotine-exposed EBs, and it has previously been shown that upregulated
   APOE leads to brain malformations and intellectual disability
   ([191]Tanaka et al., 2018). In muscle cells, increased expression of
   HMGB1 impairs cardiac excitation-contraction ([192]Zhang et al., 2014)
   and increases nicotine-induced risk for Ca^2+-associated arrhythmias
   ([193]Figure 5). In addition, the DEGs identified by scRNA-seq analysis
   are similar to bulk transcriptome studies performed in human cell lines
   ([194]Figure 2D).

   Interestingly, we found that nicotine dysregulates the cell cycle of
   endothelial, stromal, and muscle progenitor cells from G1 phase to
   S/G2M phases. Previous studies have also shown that nicotine stimulates
   the cell cycle in aortic smooth muscle cells, epithelial cells, and
   lung cancer cells ([195]He et al., 2014). One study shows that nicotine
   enhances proliferation and induces cyclin D1 to stimulate G1 to S/G2
   phase transition in human bronchial smooth muscle cells ([196]Hong
   et al., 2017), which is consistent with our findings in cluster-13
   muscle progenitor cells. A possible mechanism is that nicotine
   activates RAS/MAPK pathway via nicotinic acetylcholine receptors
   (nAChRs) to trigger a network that positively regulates cell-cycle
   progression through G1 to S, such as cyclin D ([197]Hong et al., 2017).
   Here, although no expression difference was observed on cyclin D and
   RAS/MAPK in our scRNA-seq data, we found that the components of
   cell-cycle machinery, such as HSP90AA1, TUBB4B, and TUBA1B, which are
   related to G2M transition ([198]Duggal et al., 2018), and HNRNPH1 and
   HNRNPA2B1, which are related to G1 to S transition ([199]Duggal et al.,
   2018), are upregulated in nicotine-exposed muscle cell cluster 3.

   As such, scRNA-seq analysis provides a robust tool for investigating
   cell-to-cell interactions ([200]Kawasaki et al., 2011, [201]Rehan
   et al., 2009) in development and disease pathobiology. In particular,
   an activated HMGB1-TLR4 pathway was pronounced in multiple cell types
   within EBs upon nicotine exposure. High expression of HMGB1 in multiple
   organs, perhaps induced by the secondary effects of nicotine such as
   oxidative stress, apoptosis, and inflammatory factors ([202]Loukili
   et al., 2011, [203]Scaffidi et al., 2002, [204]Kim et al., 2016), and
   is known to mediate multiple pathological conditions ([205]Ko et al.,
   2014), has been shown in smokers. For example, upregulated HMGB1
   impairs cardiac excitation-contraction coupling by enhancing
   sarcoplasmic reticulum Ca^2+ leakage through TLR4-ROS signaling in
   cardiomyocytes ([206]Zhang et al., 2014). This indicates that HMGB1
   could be a potential drug target for nicotine-induced embryonic
   defects.

   Nicotine has been suggested to mediate its function via an
   nAChR-dependent or -independent pathway. Nicotinic receptors are
   expressed in undifferentiated and differentiating cells
   ([207]Figure S3C). Studies have shown that nAChRs mediate apoptosis,
   cell proliferation, cell differentiation, regulation of intracellular
   calcium, oxidative stress, and inflammation by nicotine ([208]Dasgupta
   and Chellappan, 2006). In addition, nicotine is reported to promote
   tumor progression by binding to β-adrenergic receptors ([209]Carlisle
   et al., 2007). In addition, nicotine drives both cell proliferation and
   cell death via paracrine signaling by cell-cell interaction
   ([210]Delitto et al., 2014, [211]Scaffidi et al., 2002). Thus, nicotine
   may induce the adverse effects of EBs through nAChR-dependent and
   -independent pathways, as well as cell-cell interaction.

   Furthermore, scRNA-seq analysis can be used to optimize the treatment
   manner and period of drug use for patient-specific drug
   screening/testing. For example, with the increasing availability of
   commercial human genomic sequencing data, we can evaluate embryonic
   developmental-specific drug responses with (single nucleotide
   polymorphism) SNPs, and the response can be confirmed by gene editing
   with the CRISPR technology ([212]Seeger et al., 2017). Specifically,
   the SNP (rs141819830) of the GFI1 gene has been reported to be
   sensitive to maternal smoking in exposed neonates ([213]Gonseth et al.,
   2016). Thus, we may be able to use patient-specific hiPSC-derived EBs
   carrying this SNP to evaluate the risk of nicotine toxicity on
   embryonic development. We anticipate that this platform should help
   correlate the risk of the GFI1 SNP gene with maternal smoking during
   pregnancy.

   In summary, we used microdroplet-based scRNA-seq to investigate the
   adverse effects on heterogeneous EBs upon nicotine exposure. Our study
   offers an effective platform to evaluate the potential effects of
   nicotine on human embryonic development. Our data provide potential
   molecular mechanisms for prenatal nicotine toxicity on specific cell
   populations derived from human ESCs.

Cell Culture and Differentiation

   hESC line, H7 (WiCell Research Institute), was seeded on Matrigel (BD
   Bioscience)-coated plates in Essential 8 Medium (Thermo Fisher
   Scientific). We generated EBs by clone suspension. EBs were
   differentiated in DMEM/F12 (Gibco) supplemented with 20% FBS (Gibco),
   50 U/mL penicillin/streptomycin (Gibco), 2 mM L-glutamine (Gibco), 1×
   non-essential amino acids, and 100 μM β-mercaptoethanol (Sigma). At 90%
   confluency, hESCs were digested using 1 mL Gentle Cell Dissociation
   Reagent (STEMCELL Technologies) for 5 min. Cell clumps were pipetted
   into single cells, and 9.0 × 10^5 cells per well were seeded into
   AggreWell 800 (STEMCELL Technologies). The day after seeding,
   EBs/spheroids were harvested from AggreWell 800 plates and transferred
   into an ultra-low attachment six-well plate (Corning). Nicotine (10 μM;
   N3876, Sigma) was added into the differentiating medium, and EBs were
   fed each day for 21 days. The same volume of ethanol (459836, Sigma)
   was added into the differentiating medium as control.

scRNA-Seq Library Preparation and Analysis

   Single cells were collected from two independent EB differentiation
   experiments from day 21 EBs (control and nicotine-exposed) and
   dissociated using Accutase (STEMCELL Technologies). They were prepared
   for the single-cell library separately. In brief, cells were washed
   with 1× DPBS (Gibco) three times, strainer filtered, and re-suspended
   in 0.04% BSA. Viable single cells were loaded on to a GemCode
   Instrument (10× Genomics, Pleasanton, CA) to generate single-cell
   barcoded droplets (GEMs) using the 10× Single Cell 3′ v.2 chemistry and
   10× Chromium system as per the manufacturer's protocol. The quality of
   the resulting libraries was checked with Bioanalyzer (Agilent
   Bioanalyzer 2100). Then control and nicotine libraries were combined
   with equal molar mass and sequenced across two lanes on an Illumina
   HiSeq machine. Analyses was performed using Seurat R package. Detailed
   scRNA-seq analysis is available in the [214]Supplemental Experimental
   Procedures.

EB Diameter Measurement

   Bright-field images were captured by an SI8000 Cell Motion Imaging
   System (Sony Biotechnology) using a 4× objective. Scale was set with a
   20-cm-ruler image. The diameter of each EB was measured by drawing a
   selection line of longest distance in each EB using the Line Selection
   tool in ImageJ. Analyses were processed with ImageJ and the results
   were plotted with Prism (GraphPad).

Quantification of ROS Production and ATP

   ROS production level and ATP level in EBs were determined using
   CellTiter-Glo 2.0 (Promega) and ROS-Glo H[2]O[2] (Promega) following
   the manufacturers' instructions. In brief, EBs were exposed to nicotine
   in 96-well plates for 21 days. After treatment, H[2]O[2] substrate was
   added directly to each well of 96-well plates and incubated for 4 h.
   The resulting supernatant was collected for ROS-Glo H[2]O[2] assay and
   EBs were subsequently subjected to the CellTiter-Glo 2.0 assay to
   measure the ATP level. The luminescence intensity was measured using a
   Synergy HTX multi-mode microplate reader (BioTek).

Intracellular Calcium Imaging

   Single hESC-derived cardiomyocytes were plated onto Matrigel-coated
   coverglass (CS-24/50, Warner Instruments) at a density of ∼10,000 cells
   per square centimeter. Cells were allowed to recover for 3–4 days and
   loaded with 5 μM Fura-2 AM (Thermo Fisher Scientific) in Tyrode's
   solution (140 mM NaCl, 5.4 mM KCl, 1 mM MgCl[2], 10 mM glucose, 1.8 mM
   CaCl[2], and 10 mM HEPES; pH adjusted to 7.4 with NaOH at room
   temperature) for 30 min at room temperature. Cells were imaged on a
   customized Ti-S/L 100 Inverted Microscope-based imaging platform with a
   40× oil immersion objective (CFI Super Fluor, NA 1.30 WD 0.22). Bipolar
   pulse was used to pace cells at 0.5 and 1 HZ. Cells were kept at 37°C
   while recording. Fura-2 signals were captured in high-frame-rate video
   recording mode (512 × 512 pixels) at a speed of 50 frames per second.
   Videos were analyzed with NIS Elements Advanced Research Software
   (Nikon), and raw ratio-pair data were further processed with a
   custom-made script based on Interactive Digital Language.

Author Contributions

   H.G. designed the study, performed the experiment, analyzed the data,
   and wrote the manuscript. L.T. analyzed the scRNA-seq data and
   published bulk RNA-seq data, discussed the results, and wrote the
   manuscript. J.Z.Z. conducted Ca^2+ experiments, discussed the results,
   and revised the manuscript. T.K. performed the experiments, discussed
   the results, and revised the manuscript. D.T.P. designed the study,
   provided advice, and revised the manuscript. W.H.L. performed the
   experiments, provided advice, and revised the manuscript. J.C.W.
   designed the study, revised the manuscript, and provided funding
   support.

Acknowledgments