Cells Respond to Distinct Nanoparticle Properties ... - ACS Publications

Oct 27, 2016 - Hugh D. Mitchell,. †. Lye Meng Markillie,. †. William B. Chrisler,. † ... Jason E. McDermott,. †. Alejandro Heredia-Langner,. â...
0 downloads 0 Views 9MB Size
Download PDF
Cells Respond to Distinct Nanoparticle Properties with Multiple Strategies As Revealed by Single-Cell RNA-Seq Hugh D. Mitchell,† Lye Meng Markillie,† William B. Chrisler,† Matthew J. Gaffrey,† Dehong Hu,† Craig J. Szymanski,† Yumei Xie,† Eric S. Melby,† Alice Dohnalkova,† Ronald C. Taylor,† Eva K. Grate,† Scott K. Cooley,‡ Jason E. McDermott,† Alejandro Heredia-Langner,‡ and Galya Orr*,† †

Earth & Biological Sciences Directorate and ‡National Security Directorate, Pacific Northwest National Laboratory, Richland, Washington 99352, United States S Supporting Information *

ABSTRACT: The impact of distinct nanoparticle (NP) properties on cellular response and ultimately human health is unclear. This gap is partially due to experimental difficulties in achieving uniform NP loads in the studied cells, creating heterogeneous populations with some cells “overloaded” while other cells are loaded with few or no NPs. Yet gene expression studies have been conducted in the population as a whole, identifying generic responses, while missing unique responses due to signal averaging across many cells, each carrying different loads. Here, we applied single-cell RNA-Seq to alveolar epithelial cells carrying defined loads of aminated or carboxylated quantum dots (QDs), showing higher or lower toxicity, respectively. Interestingly, cells carrying lower loads responded with multiple strategies, mostly with up-regulated processes, which were nonetheless coherent and unique to each QD type. In contrast, cells carrying higher loads responded more uniformly, with mostly down-regulated processes that were shared across QD types. Strategies unique to aminated QDs showed strong upregulation of stress responses, coupled in some cases with regulation of cell cycle, protein synthesis, and organelle activities. In contrast, strategies unique to carboxylated QDs showed up-regulation of DNA repair and RNA activities and decreased regulation of cell division, coupled in some cases with up-regulation of stress responses and ATP-related functions. Together, our studies suggest scenarios where higher NP loads lock cells into uniform responses, mostly shutdown of cellular processes, whereas lower loads allow for unique responses to each NP type that are more diversified proactive defenses or repairs of the NP insults. KEYWORDS: single-cell sorting, transcriptional response, hierarchical clustering, functional enrichment, differential gene expression “overloaded” with hundreds or even thousands of NPs while other cells are loaded with only few or no NPs. The uneven distribution partially results from the tendency of NPs to agglomerate when the local concentration is relatively high, as is often the case in the initial exposure solution. Yet regulation of gene expression or pathways in response to NP exposure has been measured in the cell population as a whole.4 This approach identifies averaged, often most common or generic processes while leaving other critical processes undetected due to the dilution of signals across many cells, each carrying a different number of NPs. This concept was demonstrated

T

he exploding growth in the use of nanomaterials in a wide range of applications is expected to increase both the intended and unintended human exposure to engineered nanoparticles (NPs),1 but a great deal of confusion still exists about the properties that make a NP toxic or biocompatible. Distinct physical and chemical properties of the NP engage and activate distinct proteins and cellular pathways that, in turn, govern the fate of the NP and its impact on the cell and ultimately on human health. The relationships between NP properties and these key cellular processes and response are far from being understood. This knowledge gap is largely due to experimental challenges that NPs present, including the difficulty to achieve uniform NP distribution over the exposed cells. The distribution typically spans 2−3 orders of magnitude,2,3 indicating that some cells are © 2016 American Chemical Society

Received: August 12, 2016 Accepted: October 27, 2016 Published: October 27, 2016 10173

DOI: 10.1021/acsnano.6b05452 ACS Nano 2016, 10, 10173−10185

Article

www.acsnano.org

Article

ACS Nano

same QD type with different strategies, and these multiple responses, which show a coherency within a QD type, were mainly observed in cells carrying lower QD loads. Cells carrying higher loads responded more uniformly within and across QD types, with strong down-regulation of multiple functions. Interestingly, cells carrying aminated QDs at lower loads showed strong up-regulation of stress response functions, while cells carrying higher loads showed down-regulation of these functions. This observation possibly reflects a greater capacity for activating mechanisms that might help to manage adverse effects at lower QD loads versus reverting to a general shutdown of key processes at higher QD loads. In addition to a strong stress response, alternative response strategies unique to aminated QDs were observed in individual cells. One such strategy showed the same level of stress response but a marked decrease in regulatory functions over cell cycle, translation, ATP, and organelle activities. Another strategy showed an even greater stress response, coupled with up-regulations of functions involved in cell division, protein synthesis, and organelle activities. In contrast, responses unique to carboxylated QDs, observed in cells carrying lower loads, showed activation of DNA repair mechanisms, coupled with loss of regulation over mitosis, as well as increase in RNA-related activities. An alternative strategy observed in an individual cell showed an additional increase in stress responses and functions related to energy production. Considering that at very high cellular loads, aminated QDs are more toxic than carboxyated QDs, these response strategies, unique to aminated or carboxylated QDs, might be indicative of higher or lower NP toxicity, respectively. Thus, separating and focusing on individual cells carrying lower NP loads allowed us to uncover unique response patterns to distinct NP types, including multiple yet coherent response strategies. A population-based approach, however, would have obscured these response patterns within more uniform response patterns dominated by down-regulated functions in cells carrying higher NP loads.

recently using high-content imaging for linking cellular responses to intracellular QD loads at the single-cell level.5,6 Using this approach, Manshian et al. identified mechanisms in response to subtoxic QD doses, which would have been unnoticed if the behavior of the entire cell population was averaged.5,6 For example, they found that apoptosis occurs at medium to high intracellular QD loads but is impeded by autophagy at higher loads, leading to a partial recovery of cell viability. At the highest load, however, high levels of autophagy lead to cell death. Because processes that occur in a fraction of cells can be buried in the population measurement,7 unrealistically high NP doses are often used to achieve detectable responses. As a result, it is unclear whether population studies identify only generic, emergency state responses in a subset of “overloaded” cells and whether responses unique to the properties of the NPs at lower exposure and load levels, where many cells might carry not even one NP, are diluted and missed altogether. To answer these questions, we exposed alveolar epithelial cells, which present a vulnerable target for airborne NP exposure,8−10 to aminated or carboxylated quantum dots (QDs) and sorted out individual cells by their QD load for single-cell RNA-Seq analysis. QDs have been generating interest in a growing range of applications, from solar and interfacial charge transfer,11−13 including photovoltaics14 and solid-state light-emitting diodes,15 to biological research, drug delivery, and medical diagnostics.16−20 The biological and theranostic applications have gained traction with the realization that toxicity of QDs is likely to originate from the dissolution of heavy metals in the QD core,21−23 but QDs capped and coated with protective layers have the potential to be biocompatible24,25 and provide a unique tool for both treatment and imaging.19,20 However, controversy still exists about the type of coating and functional groups that could make these promising NPs biocompatible.26−30 Because of their high relevance to intended and unintended human exposure, and their bright fluorescence that allows sorting single cells by NP load, we chose to focus on aminated and carboxylated QDs to fill the gap in our understanding of the relationships between their distinct properties and the unique molecular pathways that they impact. RNA-Seq has been applied successfully to analyze the whole transcriptome in single cells, uncovering an underestimated cell-to-cell variability that was buried in cell population measurements. These studies have examined cell cycle stage assignment,31 tissue-type mapping of individual cells,32−34 conservation of transcriptome variability between species,35 stem cell detection among nonstem cells,36 reprogramming of differentiated cells to pluripotency,37 cell maturation/differentiation,38−41 detection of cancer cells among normal cells,42 gene expression variability among cells in tumors,43,44 immune response in dendritic cells,45 allele-specific expression,46 among other processes. However, none of the single-cell RNA-Seq studies thus far have attempted to identify response patterns in single cells perturbed by inherently variable external insults, such as NP exposures. Such studies have to face the challenges presented by the variability between ostensibly similar cells and, in addition, the variability introduced by the external perturbations. Here, we addressed these challenges by sorting single cells according to their NP load to minimize variability in perturbation and analyzed single-cell RNA-Seq data using multiple approaches to decipher unique response patterns to distinct NP properties. We found that cells can respond to the

RESULTS AND DISCUSSION Guided by the working hypothesis that cells loaded with a higher number of NPs enter a general emergency state, while cells loaded with a lower number of NPs respond to specific properties of the NPs, we exposed alveolar epithelial cells to a nontoxic dose of aminated or carboxylated QDs and sorted them by QD loadone cell per wellfor single-cell RNA-Seq. While available resources limited our studies to one cell type, we chose to focus on alveolar epithelial cells as they have been shown to present a vulnerable target for inhaled nanoparticles.8−10 Aminated and carboxylated cadmium selenide (CdSe) QDs with a zinc sulfide (ZnS) shell and peak emission at 655 nm were used in this study. The QD size distribution was determined from transmission electron microscopy (TEM) images (Figure 1A), where the average QD width and length were calculated by measuring 100 individual QDs for each QD type. For the aminated QDs, the average width was 9.1 ± 1.2 nm and the average length was 13.3 ± 1.4 nm. For the carboxylated QDs, the average width was 9.2 ± 1.2 nm and the average length was 13.5 ± 1.3 nm, showing no significant difference between the sizes of the two QD types. Aggregate size measured in supplemented growth medium at the concentration used in this study (10 μg/mL) showed no significant difference between the aminated and the carboxyated QDs (25.5 ± 2.3 and 26.1 ± 1.9 nm, respectively). The ζ10174

DOI: 10.1021/acsnano.6b05452 ACS Nano 2016, 10, 10173−10185

Article

ACS Nano

as well as the early time point in vivo. A nontoxic dose (10 μg/ mL) was chosen based on results using the MTS assay (Figure S1A), showing no significant decrease in cell viability, indicating no acute toxicity or extreme damage was induced in the population as a whole. It is important to note that when cells were incubated with QDs over the whole 24 h (Figure S1B), rather than 1 h, a significant cell death was detected in response to aminated QDs starting at 20 μg/mL and to carboxylated QDs starting at 50 μg/mL. This observation indicates that at very high QD loads (not pursued in our study) the aminated QDs are more toxic than the carboxylated QDs. Figure 1B shows the distribution of the relative QD load per cell (solid line), acquired by the relative fluorescence intensities of individual cells using a fluorescence activated cell sorter (FACS). The fluorescence intensity distribution spans about 2 orders of magnitude, indicating that some cells are “overloaded” with many QDs, while other cells are loaded with only few or no QDs, as shown by the overlap of the distribution with the distribution of control cells (dashed line) carrying no QDs. This distribution was used to set sorting gates (Figure 1B, vertical dotted lines) for single-cell sorting by QD load using FACS. Cells with lower QD loads were selected from the lower fluorescence intensity gate (left gate), and cells with higher QD loads were selected from the higher fluorescence intensity gate (right gate). Using confocal fluorescence microscopy, we confirmed the presence of lower and higher QD loads in the sorted cells. Figure 2 shows an example of confocal fluorescence images of alveolar epithelial cells sorted from the higher (Figure 2A) and the lower (Figure 2B) fluorescence intensity gates as well as a cell with an intensity value that is lower than the intensity in the left gate (Figure 2C). Single-Cell Transcriptomes Reveal Patterns Unique to QD Properties in Cells Carrying Lower Rather than Higher QD Loads. Individual cells, 10−13 cells for each of the five treatment groups, including lower and higher loads of aminated or carboxylated QDs and untreated control cells, were analyzed by single-cell RNA-Seq, and normalized gene expression counts were calculated for each cell. Available resources limited our study to this number of cells per treatment group, and while analyzing additional cells could increase statistical power, this number was sufficient to identify distinct transcriptional patterns for distinct treatment groups and show that certain treatments elicited more diverse, unique responses among individual cells than other treatments. Principal component analysis (PCA) was used to allow cells with similarities in gene expression profiles to cluster together. Figure 3A shows analysis using all 9233 genes that were detected in the single-cell RNA-Seq data sets. As shown in Figure 3A, PCA revealed that cells tended to cluster by QD type and load and away from control-unexposed cells (black dots). PCA also showed that cells loaded with a lower number of aminated (AL, pink) or carboxylated (CL, cyan) QDs were clustered into two distinct groups that were far apart, whereas cells loaded with a higher number of aminated (AH, red) or carboxylated (CH, blue) QDs clustered into two groups that partially overlapped with each other. These observations indicate that a lower QD load per cell evokes response patterns that are unique to the QD property, while a higher QD load evokes similar or generic response patterns regardless of QD type. PCA, in particular, visualization using the first two principal components, provides an approximation to the whole set. However, PCA can lead to incomplete conclusions regarding

Figure 1. (A) Transmission electron microscopy images of the aminated (left) and carboxylated (right) QDs that were used in this study. (B) Flow cytometry histogram of 10 000 alveolar epithelial cells exposed to a subtoxic dose of QDs (solid line). The distribution of control and unexposed cells is shown by the dashed line. The intensity of only one NP is within the background level of unexposed cells. Single cells were sorted out from the fluorescence intensity gates, marked by the vertical dotted lines. Cells with lower QD loads were taken from the left gate, and cells with higher QD loads were taken from the right gate.

potential measured in distilled water showed −0.2 ± 0.3 mV for aminated QDs and −22.83 ± 8.6 mV for carboxylated QDs. ζPotential measured in supplemented growth medium showed no significant difference between the values of the aminated (1.5 ± 3.6 mV) and carboxylated (0.44 ± 1.7 mV) QDs, possibly as the result of a protein corona that is formed on NP surfaces in supplemented growth media. This corona might be different for the carboxylated and aminated QDs, which might also contribute to the observed differences in the cellular responses. The slightly negative ζ-potential value of the aminated QDs possibly reflects a relatively low number of amine groups at the QD surface. However, the carboxylated QDs showed substantially more negative ζ-potential values compared with the values of the aminated QDs, together providing two QD types that substantially differed in the functional groups at their surface and possibly their corona but were identical otherwise. As described in Methods, we chose to expose the cells to the QDs for 1 h, after which the QDs were rinsed and the cells were incubated in growth medium for an additional 24 h before being sorted for single-cell RNA-Seq. By setting the exposure to the QDs for 1 h, the cells were loaded with a discrete QD load per cell, rather than the continuous increase in QD load that would have occurred over longer exposure periods, with cellular processes also shifting over time. The 24 h response time was selected as it is commonly used in NP exposure in vitro studies, 10175

DOI: 10.1021/acsnano.6b05452 ACS Nano 2016, 10, 10173−10185

Article

ACS Nano

Figure 2. Confocal fluorescence images in 3D of cells taken from the higher (A) and lower (B) intensity gates (right and left gates in Figure 1B, respectively), as well as a cell with a lower intensity value than the values in the left gate (C). Membrane-specific fluorescent probe was used to outline the membrane (green). By removing or dimming the green color in images 2 and 3, respectively, the QDs (red or yellow) inside the cells can be seen, with few QDs found at the cell surface (examples indicated by the arrows). The particles seem elongated because the lateral resolution is twice as good as the axial resolution in confocal fluorescence microscopy. Scale bars: 4 μm.

than all 9233 genes found in the RNA-Seq data sets. This separation is demonstrated by the PCA shown in Figure 3B. Conventional differential expression analysis allowed us to identify genes that were expressed at statistically different levels in cells from each treatment group compared to control cells. Differentially expressed genes were then submitted to functional enrichment analysis to determine functional pathways and categories likely to be over-represented among differentially expressed genes. It should be noted that increases in gene expression for a subset of genes involved in a given pathway/function do not necessarily translate to up-regulation of this pathway/function because changes in pathway antagonists will have opposite effects. The enrichment results using all 9233 detected genes were rolled-up to a functional cluster level (see Methods) and represented in a heatmap according to a significance score for each cluster (Figure 4). As Figure 4 demonstrates, the two higher QD load groups (AH and CH) mostly showed downregulated responses, while the lower aminated QD load (AL) group showed the greatest number of up-regulated responses. Some of these (heat shock protein, lysosome, endoplasmic reticulum (ER), unfolded protein binding, apoptosis, and cell redox homeostasis) were uniquely up-regulated by the AL

which variables are most relevant to the distribution of the data points in the two-dimensional space when the first two principal components are used. This can be a result of very high within-class variability, highly correlated variables, or a combination of the two.47 One way to address this problem when dealing with large data sets is to systematically select a subset of variables that adequately represent the entire data set. Classification trees can be useful in this context because these algorithms select, using a greedy approach, only those variables that can best discriminate between groups or classes while ignoring the rest. A single classification tree will, usually, select only a few variables useful for discrimination. Other variables can be selected by applying classification trees to the data in a recursive way. To identify the functions or processes that can separate most clearly the responses of all treatment groups, we therefore used classification trees that optimize classifications of the individual cells into their respective treatments. Using this approach, we identified the top 100 genes that contributed the most to the separation between the treatment groups. We found that these top 100 genes (listed in Table S1), chosen by the classifier, could more effectively segregate cells by their treatment groups 10176

DOI: 10.1021/acsnano.6b05452 ACS Nano 2016, 10, 10173−10185

Article

ACS Nano

aminated QDs (Figure S1B), these unique responses by the AL group, dominated by up-regulations of stress responses and certain fine-tuning functions, might reflect distinct response patterns to more toxic perturbations. Functional analysis of the top 100 genes, selected by classification trees as described above, revealed regulations of two main functional categories, including cell cycle and stress response. Genes relating to these two functional categories from the identified 100 are listed in Table S2. Pathways connected to these two categories were highly prevalent in the functional analysis described in Figure 4 using all 9923 genes. These functions, indicated with red stars in Figure 4, represent pathways with unique signatures across the four treatment conditions, consistent with expectations that different QD types and loads will cause cells to modulate their growth and stress response differently, depending on the severity of the insult. For example, the CL treatment group was differentiated from all other treatment groups by showing mixed up- and downregulation of cell-cycle-related genes and normal levels of chromosome kinetochore functions, whereas the other three treatment groups showed clear up-regulation of these functions (Figure 4, red stars). Further, the CL group uniquely showed strong up-regulation of the microtubule organizing center and DNA replication functions (Figure 4, red stars), where the former may represent perturbation of key regulatory components that repress mitosis through perturbation of cytoskeletal machinery, and inspection of genes from the latter reveals a strong component of DNA repair genes. Other functions uniquely up-regulated by the CL group included RNA binding and transport. Considering the lower toxicity of the carboxylated QDs (Figure S1B), these unique responses by the CL group, including activation of DNA repair mechanisms, coupled with loss of regulation over mitosis as well as increase in RNA-related activities, might reflect distinct response patterns to perturbations with lower toxicity levels. Individual Cells Can Respond with Different Strategies That Together Create a Unique Pattern of Response to a Distinct QD Property. To identify response patterns or strategies within each treatment group, we applied hierarchical clustering to all 9233 detected genes and plotted the various gene expression responses with high variance across treatments and cells (Figure 5). The hierarchical clustering in Figure 5A shows that many genes exhibit extreme (high or low) expression ranges in a subset of cells. These divergent cells (indicated with red bars in Figure 5A) were present only in the lower QD load or control groups and not in the higher QD load groups. This observation is in line with our finding that the responses of cells loaded with a higher number of QDs tend to be more similar, regardless of QD type, while responses of cells loaded with a lower number of QDs tend to be more unique to QD type. One possible result of radically diverse responses is loss of treatment group cohesiveness, such that some divergent cells could fail to bear any resemblance to other cells within their group. Interestingly, however, hierarchical clustering of these divergent cells alone (Figure 5B) reveals that they still cluster according to treatment, as shown by the dendrogram on the left, indicating a degree of coherency in their divergent behavior. To understand the divergent behavior of cells within their treatment group in a quantitative way, we applied hierarchical clustering for each group separately (Figure 6). Figure 6 shows the hierarchical clustering height of each treatment group, where the emerging patterns reveal that the lower QD load

Figure 3. Principal component analysis of the single-cell RNA-Seq data, using normalized gene expression counts for each cell. Each single-cell data for each treatment group is numbered and color coded. (A) First two principal components (x- and y-axes) for the normalized gene expression counts of all 9233 detected genes. (B) First two principal components for the normalized gene expression counts of the top 100 genes, selected using classification trees by their effectiveness to separate treatment groups. Key: Controlunexposed cells (black); AL-aminated QDs at a lower load (pink); CL-carboxylated QDs at a lower load (cyan); AH-aminated QDs at a higher load (red); CH-carboxyated QDs at a higher load (blue).

group while showing no or down-regulation by the other groups, indicating a strong stress response unique to the AL group. In contrast, several of these and other stress-related processes were down-regulated in the higher aminated QD load (AH) group. This result is surprising because AL cells carried lower QD loads than AH cells, yet they appear to be responding with stress pathways to a greater extent than AH cells. This difference might reflect a responsive cellular state with a potential for repair and survival of the AL cells versus an “overwhelmed” cellular state of the AH cells, where certain cellular machineries might have been shut down. Another function uniquely up-regulated by the AL group was the ubl conjugation function. This function represents ubiquitinationrelated genes involved in cell cycle checkpoints (Anapc11, Bard1, Bub1b, Ccnb1, Ccnd1, Cksn1a, Cdkn1b, Trim37, Uhrf2), proliferation (Birc5, Ctnnb1, Lig1, Nusap1, Sgol1), and chaperone activity (Cct8, Tmbim6). Activation of these pathways in the AL cells shows an emphasis on the fine-tuning machinery, which was missing in cells carrying higher loads of these same QDs. Considering the higher toxicity of the 10177

DOI: 10.1021/acsnano.6b05452 ACS Nano 2016, 10, 10173−10185

Article

ACS Nano

Figure 4. Group-wise functional enrichment. Up- and down-regulated genes from differential expression analysis using all 9233 detected genes were submitted to functional enrichment analysis separately to determine the direction of change in functions/pathways. Enrichment results were clustered, matched across treatment groups, and assigned scores. Functional enrichment of down-regulated genes were multiplied by −1. Heatmap units with two colors indicate that both up- and down-regulated genes were found to be enriched in a given function/pathway. Stars indicate functions/pathways related to cell cycle and stress response identified using the top 100 genes, selected by their effectiveness to separate treatment groups (Tables S1 and S2 and Figure 3B). Treatment group abbreviations are the same as that in Figure 3.

are listed in Tables S3−S6. Figure S2 shows cluster dendrograms using the identified genes, and Figure 7 shows a functional heatmap of the pathways/functional categories found to be enriched in the identified 200 gene lists. Interestingly, cell 7 in the AL treatment group showed an even greater stress response than other cells in this group, while also up-regulating other functions, including actin and microtubule processes involved in cell division and organelle activity, as well as in protein synthesis (translation and polymerization) and ubiquitination, RNA binding, and ATP synthesis. Cell 8 in the AL treatment group took yet another strategy, showing stress responses at the same level as the rest of the cells in the group but a marked decrease in the bilateral symmetry function, which includes genes likely to oppose or regulate cell cycle (Atm, Fkbp8, Notch2), as well as in processes related to ATP binding, translation initiation, and organelle activity. Cell 10 in the CL group showed increased stress responses compared with other cells in its treatment group, with even higher activities in unfolded protein binding, cell redox homeostasis, and oxidative stress response, as well as in functions related to energy production. Cell 2 in the CH group showed marked down-regulation in transcription-related

groups (AL and CL) are highly variable, while the higher load of aminated QD (AH) group shows a markedly lower degree of variation. Since technical noise must be present to similar degrees in all treatment groups, we reasoned that the variability in the treatment group with the least degree of variation (AH) represents the maximum possible effect of technical noise, and that additional variations stem from true biological differences. To identify and characterize the divergent responses within treatment groups, we therefore chose a cutoff well above the maximum clustering height of the AH treatment group (Figure 6, red line). Cells with branch points above this cutoff were considered to have significantly different behaviors from the rest of their group. Using this cutoff, the AL, CL, and CH treatment groups were subdivided into two groups, above and below the cutoff value, and the 200 genes most responsible for the segregation were identified. This was done by searching the 9233 identified genes for those genes that produced maximal separation between the target cell (above the cutoff line) and the cells that show common behavior within the treatment group, scaled by the standard deviation of the gene expression counts in the common group. Top genes identified for each treatment group 10178

DOI: 10.1021/acsnano.6b05452 ACS Nano 2016, 10, 10173−10185

Article

ACS Nano

capacity to address the perturbation in the divergent cells, compared with the rest of the cells in their group. Two additional questions come to mind. First, would an intermediate QD load elicit a distinct transcriptional response, or would it elicit a response that is similar to the responses elicited by either the higher or lower QD loads, suggesting a response threshold? Second, would narrower sorting gates minimize the variability that we found between cells within a treatment group? To address these questions, we exposed the cells to carboxylated QDs at a nontoxic dose (10 μg/mL) and sorted single cells from three distinct gates (lower, intermediate, and higher loads) that were nearly 3-fold narrower than the gates used earlier (Figure S3A). For each QD load group, as well as control, 11 cells were analyzed by single-cell RNA-Seq. These carboxlated QDs had a slightly more negative ζ-potential value (−24.5 ± 6.2 mV) and identical core/shell structure and size distribution as the carboxylated QDs used earlier. Figure 8A shows the PCA plot of the transcriptional response to the three QD load levels compared to control, using normalized gene expression values for the 10 457 genes that were detected in all treatment groups. The clustering pattern shows that each QD load level elicited a distinct response, each of which is distinct from control cells. This is further demonstrated in Figure 8B, which depicts the overlap in differentially expressed genes for each QD load group (compared to control). Each treatment group clearly expressed a large proportion of genes that were distinct from the genes expressed by the other two groups. Functional enrichment analysis revealed responses unique to cells carrying an intermediate QD load, including those related to Golgi and ER components, nucleotide binding, protein folding, and apoptosis. These observations suggest that each QD load level elicits unique responses and argue against the presence of a transcriptional response tipping point that determines which of two responses a system will display. To determine whether the narrower sorting gates led to lower variability within QD load groups, we performed hierarchical clustering within each treatment group as before. Figure 8C shows relatively tight clustering in control cells and cells carrying higher QD loads. However, the intermediate and lower load groups displayed looser clustering, thus more widely varying behavior among cells. These results replicate the trends observed using wider sorting gates (Figure 6), showing that cells carrying lower QD loads respond with multiple strategies while cells carrying higher loads respond more uniformly. Importantly, these results demonstrate that lower deviation in QD load still allowed for biological variability within treatment groups. Since less variation between cells within a treatment group is expected to cause lower standard deviations for individual genes, we compared gene variance in treated cells, sorted using wider or narrower gates, to their respective controls. Figure S3B shows median standard deviation ratios of treatment versus control groups and reveals that narrower sorting gates did not strikingly lower gene expression variability. Thus, the alternative response strategies that we found in individual cells within a treatment group are likely due to predominantly biological variability rather than variation in QD load.

Figure 5. Hierarchical clustering analysis of the single-cell RNA-Seq data, using normalized gene expression counts for each cell. (A) All 9233 genes are represented in the heatmap, where the color codes indicate normalized count levels. Rows are individual cells, while columns are individual genes. Genes underwent hierarchical clustering, while rows are grouped according to treatment. Red bars delineate cells with divergent gene expression behavior compared to other cells in their treatment group. (B) All 9233 genes are represented for the nine divergent cells (red bars in A), with all other parameters the same as in A. Hierarchical clustering analysis shows that these cells still cluster within their treatment groups, as indicated by the dendrogram on the left. Treatment group abbreviations are the same as that in Figure 3.

processes, indicating a more profound shutdown of certain cellular processes compared with the rest of the cells in its higher QD load group. Thus, individual cells within a treatment group can manifest significantly distinct behaviors which can be completely missed due to signal averaging that occurs in population-based studies. The contribution of a minority of divergent cells to the population is unknown but may have significant impact through paracrine effects as has been shown previously.45 Interestingly, the only shared process across all divergent cells was the upregulation of cellular respiration (electron transport; mitochondrial membrane) compared with their treatment group (Figure 7). This up-regulated function, which was down-regulated in the higher load groups of both QD types, might reflect a greater

CONCLUSION The main findings to emerge from this work are that individual cells within a population can respond with different biological 10179

DOI: 10.1021/acsnano.6b05452 ACS Nano 2016, 10, 10173−10185

Article

ACS Nano

Figure 6. Dendrograms of the hierarchical clusters obtained using normalized gene expression counts for all 9233 genes for all treatment groups. The y-axis shows the distance between connected points. The horizontal red line represents the cutoff line, where cells with branching points above this value were considered to have significantly different behaviors compared with the rest of the cells in their group. Cell numbers and treatment group abbreviations are the same as that in Figure 3.

processes to the exact same QD type and load, and these multiple yet coherent response strategies create a unique pattern of response to a distinct QD property. Furthermore, these unique responses to distinct QD properties are revealed mainly in cells that are loaded with a lower number of QDs, while cells carrying higher loads respond more uniformly within and across the two QD types, with a marked down-regulation of multiple processes. Counterintuitively, we found that cells carrying lower loads of aminated QDs showed strong upregulation of stress responses, which were markedly downregulated in cells carrying higher loads of the same aminated QDs. This observation suggests that cells carrying lower loads may have greater capacity to proactively manage the adverse effects and elicit defense responses, such as heat shock and oxidative processes, which can effectively address the unique challenges presented by the QDs. In contrast, cells carrying higher QD loads may be “resigned” to an unfavorable outcome and therefore manage the extreme perturbation mainly with a general shutdown of key processes. Higher loads of either aminated or carboxylated QDs led to a general shutdown of multiple processes, including translation, stress response, cellular transport, and energy production, indicating what appears to be a general failure to appropriately respond to the perturbation. Interestingly, we also found that cells carrying an intermediate QD load elicited a unique response in comparison

with higher and lower loads, ruling out a threshold or a tipping point in the transcriptional response with increasing load. By recursive applications of classification trees to the singlecell RNA-Seq data, we found that the genes that can most effectively distinguish between the treatment groups were highly enriched in stress responses and cell cycle processes. This observation indicates that cells responded to each QD type and load by uniquely regulating these two functional categories. Considering that at high cellular loads, aminated QDs are more toxic than carboxyated QDs, it is possible that upregulation of stress responses, with or without additional regulation of functions related to cell division, protein synthesis, and organelle activities, seen clearly in cells carrying lower loads of aminated QD, is indicative of a cellular attempt to recover from toxic NP properties. In contrast, up-regulation of DNA repair mechanisms and RNA activities, coupled with decreased regulation of mitosis, with or without additional increase in stress responses and energy production functions, which were seen clearly in cells carrying lower loads of carboxylated QDs, might be indicative of a viable response strategy to NPs with lower toxicity levels. Thus, separating and focusing on individual cells carrying lower NP loads allowed us to uncover unique response patterns to distinct NP properties, including multiple yet coherent response strategies. Multiple response strategies within the lower QD load group were also preserved 10180

DOI: 10.1021/acsnano.6b05452 ACS Nano 2016, 10, 10173−10185

Article

ACS Nano

also shifting over time. A nontoxic dose (10 μg/mL) was chosen based on results using the MTS assay as described below. Toxicity Assessment. Cell viability was determined using the MTS proliferation assay (CellTiter 96, Promega) as described previously.50,51 Briefly, measurements were performed on three separate days in four 96-well plates, with a total of 12 wells per QD concentration. For each experiment, 3000 cells were seeded in each well for 48 h incubation at 37 °C. QDs were then added at the indicated concentrations for 1 h incubation at 37 °C, after which the cells were washed and incubated for additional 24 h. The cells were then assayed following the manufacturer’s protocol. Cells that were incubated with the QDs for the whole 24 h were also assayed to allow results to be compared with other studies applying this commonly used incubation protocol. Cell viability values were determined by dividing each well absorbance value by the averaged absorbance of positive control (lysed cells following the manufacturer protocol, no MTS reaction) and subtracting from 1 to convert from fraction of dead cells to fraction of live cells. These individual normalized values were then averaged for each QD concentration, and a standard deviation was determined. Significance was calculated by two-tailed Student’s t test. Nanoparticle Characterization and Preparation of QD Suspensions. Amine- and carboxyl-functionalized QDs with peak emission at 655 nm were purchased from Molecular Probes (Life Technologies). These QDs have a core made of cadmium selenide (CdSe), a semiconductor shell made of zinc sulfide (ZnS) to stabilize the core, and an amphiphilic polymer coating to which carboxyl groups are directly coupled or amine groups are coupled via polyethylene glycol linkers. The QD size distribution was determined from images acquired by TEM. The QD stock solution was diluted to 1:50 in dH2O, and 5 μL aliquots were applied onto a 100 Cu mesh grid coated with Formvar and carbon (Electron Microscopy Sciences), gently whisked off with an edge of a filter paper. The grid was allowed to airdry, and samples were examined with the Tecnai T-12 TEM (FEI) operating at 120 kV with a LaB6 filament. Images were collected digitally with a 2 × 2k Ultrascan 1000 CCD (Gatan). Aggregate size distribution was measured in supplemented growth medium (RPMI supplemented with FBS) using dynamic light scattering (ZetaPALS, Brookhaven Instruments), and ζ-potential was measured in distilled water and supplemented growth medium using a zeta potentiometer (ZetaPALS, Brookhaven Instruments). QDs were suspended in complete cell culture medium (RPMI, 10% FBS, 0.1% penicillin/streptomycin) at the indicated concentrations. These solutions were then sonicated with a cup horn sonicator (Branson Digital Sonifer 450) for 2 × 30 s with a 30 s delay between with 400 W on ice. The solutions were then warmed in a 37 °C bath for 10 min before introduction onto cells. QD suspensions were prepared immediately before exposure. Single-Cell Analysis and Sorting by Fluorescence Intensity. Following 1 h incubation with the QDs and an additional 24 h incubation with no QDs, the cells were dissociated off the plate by incubating with 0.25% trypsin in EDTA (Life Technologies) in growth medium for 5 min followed by two washes in PBS. This was done by centrifugation at 1000 rpm for 3 min, discarding the supernatant and resuspending the cells. The fluorescence activated cell sorter (Influx, BD Biosciences) was used for both analysis and sorting. Excitation was done using 355 nm, and emission was acquired at 670/30 nm wavelengths. The fluorescence intensity of individual cells (10 000 cells per QD type) were measured and plotted to create the population distribution and identify the sorting gates for the lower and higher QD loads. Once sorting gates were identified, single cells were sorted into a 96-well plateone cell per wellcontaining RNAlater (Life Technologies). Fluorescence Imaging. Confocal fluorescence microscopy (LSM 710, Zeiss) was used to confirm the lower and higher QD loads in selected cells. The cells were stained with 50 μg/mL WGA-Alexa 488 (Life Technologies) for 10 min to outline the cell membrane. Both red-emitting QDs (Qdot 655) and green fluorescence WGA-Alexa 488 were excited by 488 nm wavelength laser simultaneously. Their fluorescence emission was acquired separately into two channels,

Figure 7. Functional enrichment in divergent cells. Functional enrichment analysis was conducted on the top 200 genes found to strongly differentiate divergent cells from other cells in their treatment group (listed in Tables S3−S6). Heatmap values are colored according to the significance of the ontology category enrichment (color shade) and by the direction of change compared to the remaining cells in the group (blue = down-regulation, red/ orange = up-regulation). Cell numbers and treatment group abbreviations are the same as that in Figure 3.

when using narrower sorting gates, supporting the idea that these differences are likely due to biological variability, rather than variation in QD loads. These response patterns would have been buried otherwise within more uniform response patterns dominated by down-regulated functions in cells carrying higher NP loads, had a traditional population-based approach been used. In summary, our studies suggest a scenario where higher NP loads lock cells into a somewhat uniform response involving a marked shutdown of cellular processes, whereas lower loads allow for responses that are unique to the NP types and are more diversified, proactive defenses against the NP insult. Future single-cell studies are likely to identify additional alternative response strategies to the same NP type, allowing for a more complete comprehension of the complex relationships between NP properties and cellular response.

METHODS Cell Culture and QD Exposure. A nontumorigenic alveolar type II epithelial cell line (C10), derived from a healthy lung of an adult mouse, was used in this study. The cells were grown and exposed to the NPs as described earlier.48,49 Briefly, the cells were grown in RPMI growth medium supplemented with 10% FBS, 2 mM L-glutamine, 100 U/mL penicillin, and 100 μg/mL streptomycin. A total of 3 × 105 cells were seeded on cell culture-treated 35 mm plastic dishes (BD Falcon, 11.78 cm2 effective growth area). The cells were kept in culture at 37 °C for about 3 days until full confluence was reached, after which the cells were exposed to the QDs at 10 μg/mL suspended in 1.3 mL growth medium for 1 h at 37 °C. Cells were then rinsed and incubated in growth medium for an additional 24 h before sorting for single-cell RNA-Seq. The 1 h incubation with the QDs allowed a more discrete QD load per cell to be achieved, in contrast to the 24 h incubation where the load would continuously increase, with cellular processes 10181

DOI: 10.1021/acsnano.6b05452 ACS Nano 2016, 10, 10173−10185

Article

ACS Nano

Figure 8. Comparison of single-cell transcriptional responses across low, intermediate, and high QD loads, in cells selected using narrower sorting gates for each load group. (A) PCA plot using normalized gene expression counts of all 10 457 genes detected by single-cell RNA-Seq for each cell (numbers), showing distinct transcriptional responses for each of the three QD load groups. (B) Overlap of differentially expressed genes for each QD load group compared to control, showing that each group expressed a large proportion of genes that were distinct from the genes expressed by the other two groups. (C) Cluster dendrograms showing degree of variability (distance) between single cells within each QD load group, replicating the trends observed using wider sorting gates, where cells carrying lower QD loads responded with multiple strategies while cells carrying higher loads responded more uniformly. where the QD emission was collected at 636−703 nm, and the Alexa 488 emission was collected at 498−586 nm. Image analysis and 3D visualization were done using Volocity software (PerkinElmer). Single-Cell RNA-Seq. cDNA libraries from single cells were generated using the Smart-Seq protocol,50 which takes advantage of a reverse transcriptase from the Moloney murine leukemia virus to generate and amplify full-length cDNA from single cells. This method enriches for transcripts with an intact 5′ end and eliminates the need for second strand synthesis. This approach has been shown to expand coverage along the entire length of mRNA transcripts with starting mRNA as low as 10 pg.52−54 The Smart-Seq kit (SMARTer ultra low input RNA kit, Clontech) was used to construct full-length cDNA according to the manufacturer’s instruction. The full-length cDNA was amplified and fragmented using NEBNext dsDNA Fragmentase (New England Biolabs) followed by ligation with the SOLiD adapters

according to the manufacturer’s recommended protocol. cDNA fragments, 150−250 bp in size, were size-selected using Agencort Ampure XP (Beckman Coulter). The isolated cDNA was amplified through 15 amplification cycles to produce enough templates for emulsion clonal bead amplification using the SOLiD EZ Bead system (Life Technologies). The 50 base read sequences produced by the 5500xl SOLiD platform were mapped in color space using LifeScope software onto the mouse reference genome (UCSC Genome Bioinformatics http://hgdownload.cse.ucsc.edu/downloads. html#mouse). The software generates a gene count file, where the sum of the base counts across the entire gene length and RPKM values are provided for each gene, and BAM files containing the sequence of all mapped reads and their location are generated. A summary report for genomic alignment rates and filtering statistics is also generated for evaluating the quality of the libraries. Data were examined for rRNA 10182

DOI: 10.1021/acsnano.6b05452 ACS Nano 2016, 10, 10173−10185

Article

ACS Nano contamination and 3′ or 5′ coverage bias, and negative control libraries (empty wells) were included. Bioinformatics and Statistical Approaches Applied to the Single-Cell RNA-Seq Data. Hierarchical clustering, PCA, heatmap generation, and differential expression analysis were carried out using the R statistical programming platform, with differential expression analysis performed using the biocondoctor package edgeR.55 Functional enrichment results were acquired using the DAVID web resource56 to identify the most significantly regulated functions or pathways affected by NP exposure. The DAVID functional annotation tool utilizes the Fisher Exact test to measure enrichment in biological functions for all significantly regulated transcripts. The “functional annotation clustering” option was selected for DAVID output, which results in significant hits being grouped according to gene hit overlap between enrichment categories. To obtain the summarized enrichment results in Figures 4 and 7, enrichment hit clusters were matched across samples when 50% of the category names within a cluster were identical across samples. Hit clusters were automatically given names based on the most common words in the enrichment hit names. Some names were manually amended for clarity. Up- and down-regulated gene lists were submitted separately, so that enrichment of both upand down-regulated activity could be assessed. Enrichment scores (assigned by DAVID) from each cluster using down-regulated genes were multiplied by −1 to allow regulation in both directions to be displayed in the heatmap. For variance analysis (Figure S3B), genespecific standard variations were collected across single cells for each condition (i.e., resulting in a standard deviation value for each gene and each treatment group). Bar plots in Figure S3B represent the ratio of (1) the mean of all genes’ standard deviations for the treatment group to (2) the mean of all genes’ standard deviations for the control. Error bars are 95% confidence intervals from a bootstrap approach, in which the above ratio was calculated 10 000 times by sampling control and treatment gene-specific standard deviations with replacement. Confidence intervals are the middle 95% of means from the bootstrap procedure. Classification Trees. To identify genes most responsible for separating individual cells from the rest of the group, classification trees were applied recursively to the set of gene expression data, removing each time variable already selected from consideration. This approach enables searching for the top 100 genes that best discriminate between the five groups of cells. Functional Discrimination of Individual Cells. To identify genes most responsible for separating individual cells from the rest of the group, each gene was scored according to the following formula, where targ = the expression level of the target cell and nontarg = the mean expression level of all other cells in the treatment group.

rank =

ACKNOWLEDGMENTS This work was supported by PNNL’s Signature Discovery Initiative, and by the National Institute of Environmental Health Sciences (1RC2ES018786-01 to G.O). The research was performed using the Environmental Molecular Sciences Laboratory (EMSL), a national scientific user facility sponsored by the Department of Energy’s Office of Biological and Environmental Research and located at PNNL. REFERENCES (1) Roco, M., Bainbridge, W. S., Eds. Nanotechnology: Societal Implications IMaximizing Benefits for Humanity; Springer: Dordrecht, The Netherlands, 2007. (2) Orr, G.; Chrisler, W. B.; Cassens, K. J.; Tan, R.; Tarasevich, B. J.; Markillie, L. M.; Zangar, R. C.; Thrall, B. D. Cellular Recognition and Trafficking of Amorphous Silica Nanoparticles by Macrophage Scavenger Receptor A. Nanotoxicology 2011, 5, 296−311. (3) Orr, G.; Panther, D. J.; Cassens, J.; Phillips, J. L.; Tarasevich, B. J.; Pounds, J. G. Syndecan-1 Mediates the Coupling of Positively Charged Submicrometer Amorphous Silica Particles with Actin Filaments Across the Alveolar Epithelial Cell Membrane. Toxicol. Appl. Pharmacol. 2009, 236, 210−220. (4) Tilton, S. C.; Karin, N.; Tolic, A.; Xie, Y.; Lai, X.; Hamilton, R. F.; Waters, K.; Holian, A.; Witzmann, F. A.; Orr, G. Three Human Cell Types Respond to Multi-walled Carbon Nanotubes and Titanium Dioxide Nanobelts with Cell-specific Transcriptomic and Proteomic Expression Patterns. Nanotoxicology 2014, 8, 533−548. (5) Manshian, B. B.; Munck, S.; Agostinis, P.; Himmelreich, U.; Soenen, S. J. High Content Analysis at Single Cell Level Identifies Different Cellular Responses Dependent on Nanomaterial Concentrations. Sci. Rep. 2015, 5, 13890−13897. (6) Manshian, B. B.; Abdelmonem, A. M.; Kantner, K.; Pelaz, B.; Klapper, M.; Nardi Tironi, C.; Parak, W. J.; Himmelreich, U.; Soenen, S. J. Evaluation of Quantum Dot Cytotoxicity: Interpretation of Nanoparticle Concentrations Versus Intracellular Nanoparticle Numbers. Nanotoxicology 2016, 10, 1318−1328. (7) Wang, D.; Bodovitz, S. Single Cell Analysis: The New Frontier in ‘Omics.’. Trends Biotechnol. 2010, 28, 281−290. (8) Oberdorster, G.; Oberdorster, E.; Oberdorster, J. Nanotoxicology: An Emerging Discipline Evolving from Studies of Ultrafine Particles. Environ. Health Perspect. 2005, 113, 823−839. (9) Donaldson, K.; Borm, P. J. A.; Oberdorster, G.; Pinkerton, K. E.; Stone, V.; Tran, C. L. Concordance Between in vitro and in vivo Dosimetry in the Proinflammatory Effects of Low-toxicity, Lowsolubility Particles: The Key Role of the Proximal Alveolar Region. Inhalation Toxicol. 2008, 20, 53−62. (10) Mercer, R. R.; Hubbs, A. F.; Scabilloni, J. F.; Wang, L. Y.; Battelli, L. A.; Schwegler-Berry, D.; Castranova, V.; Porter, D. W. Distribution and Persistence of Pleural Penetrations by Multi-walled Carbon Nanotubes. Part. Fibre Toxicol. 2010, 7, 28−36. (11) Tvrdy, K.; Frantsuzov, P. A.; Kamat, P. V. Photoinduced Electron Transfer from Semiconductor Quantum Dots to Metal Oxide Nanoparticles. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 29−34. (12) Zhu, G.; Pan, L. K.; Xu, T.; Sun, Z. One-Step Synthesis of CdS Sensitized TiO(2) Photoanodes for Quantum Dot-Sensitized Solar Cells by Microwave Assisted Chemical Bath Deposition Method. ACS Appl. Mater. Interfaces 2011, 3, 1472−1478. (13) Zhu, G.; Pan, L.; Xu, T.; Zhao, Q. F.; Lu, B.; Sun, Z. Microwave Assisted CdSe Quantum Dot Deposition on TiO(2) Films for Dyesensitized Solar Cells. Nanoscale 2011, 3, 2188−2193. (14) Nimmo, M. T.; Caillard, L. M.; De Benedetti, W.; Nguyen, H. M.; Seitz, O.; Gartstein, Y. N.; Chabal, Y. J.; Malko, A. V. Visible to Near-Infrared Sensitization of Silicon Substrates via Energy Transfer from Proximal Nanocrystals: Further Insights for Hybrid Photovoltaics. ACS Nano 2013, 7, 3236−3245. (15) Kim, S.; Im, S. H.; Kim, S. W. Performance of Light-emittingdiode Based on Quantum Dots. Nanoscale 2013, 5, 5205−5214.

targ − mean(nontarg) st dev(nontarg)

Genes were ranked according to these scores, and the top genes were selected for functional enrichment analysis.

ASSOCIATED CONTENT S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acsnano.6b05452. Supplementary Figures S1−S3 and Tables S1−S6 (PDF)

AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest. 10183

DOI: 10.1021/acsnano.6b05452 ACS Nano 2016, 10, 10173−10185

Article

ACS Nano (16) Zhao, J.; Qiu, X.; Wang, Z.; Pan, J.; Chen, J.; Han, J. Application of Quantum Dots as Vectors in Targeted Survivin Gene siRNA Delivery. OncoTargets Ther. 2013, 6, 303−309. (17) Knipe, J. M.; Peters, J. T.; Peppas, N. A. Theranostic Agents for Intracellular Gene Delivery with Spatiotemporal Imaging. Nano Today 2013, 8, 21−38. (18) Kim, S.; Lim, Y. T.; Soltesz, E. G.; De Grand, A. M.; Lee, J.; Nakayama, A.; Parker, J. A.; Mihaljevic, T.; Laurence, R. G.; Dor, D. M.; Cohn, L. H.; Bawendi, M. G.; Frangioni, J. V. Near-infrared Fluorescent Type II Quantum Dots for Sentinel Lymph Node Mapping. Nat. Biotechnol. 2004, 22, 93−97. (19) Chong, L.; Vannoy, C. H.; Noor, M. O.; Krull, U. J. Intracellular Nucleic Acid Interactions Facilitated by Quantum Dots: Conceptualizing Theranostics. Ther. Delivery 2012, 3, 479−499. (20) Wu, Y.; Li, L.; Mao, Y.; Lee, L. J. Static Micromixer-coaxial Electrospray Synthesis of Theranostic Lipoplexes. ACS Nano 2012, 6, 2245−2252. (21) Chen, N.; He, Y.; Su, Y.; Li, X.; Huang, Q.; Wang, H.; Zhang, X.; Tai, R.; Fan, C. The Cytotoxicity of Cadmium-based Quantum Dots. Biomaterials 2012, 33, 1238−1244. (22) Zaitseva, N.; Manna, L.; Gerion, D.; Saw, C. K. Precipitation of Selenium from CdSe Nanocrystal Solutions. Adv. Mater. 2005, 17, 1321−1324. (23) Rikans, L. E.; Yamano, T. M. Mechanisms of Cadmiummediated Acute Hepatotoxicity. J. Biochem. Mol. Toxicol. 2000, 14, 110−117. (24) Kirchner, C.; Liedl, T.; Kudera, S.; Pellegrino, T.; Muñoz Javier, A.; Gaub, H. E.; Stölzle, S.; Fertig, N.; Parak, W. J. Cytotoxicity of Colloidal CdSe and CdSe/ZnS Nanoparticles. Nano Lett. 2005, 5, 331−338. (25) Hoshino, A.; Fujioka, K.; Oku, T.; Suga, M.; Sasaki, Y. F.; Ohta, T.; Yasuhara, M.; Suzuki, K.; Yamamoto, K. Physicochemical Properties and Cellular Toxicity of Nanocrystal Quantum Dots Depend on Their Surface Modification. Nano Lett. 2004, 4, 2163− 2169. (26) Stern, S. T.; Zolnik, B. S.; McLeland, C. B.; Clogston, J.; Zheng, J.; McNeil, S. E. Induction of Autophagy in Porcine Kidney Cells by Quantum Dots: A Common Cellular Response to Nanomaterials? Toxicol. Sci. 2008, 106, 140−152. (27) Ho, C. C.; Chang, H.; Tsai, H. T.; Tsai, M. H.; Yang, C. S.; Ling, Y. C.; Lin, P. Quantum Dot 705, a Cadmium-based Nanoparticle, Induces Persistent Inflammation and Granuloma Formation in the Mouse Lung. Nanotoxicology 2013, 7, 105−115. (28) Lin, P.; Chen, J. W.; Chang, L. W.; Wu, J. P.; Redding, L.; Chang, H. Computational and Ultrastructural Toxicology of a Nanoparticle, Quantum Dot 705, in Mice. Environ. Sci. Technol. 2008, 42, 6264−6270. (29) Zhang, T.; Stilwell, J. L.; Gerion, D.; Ding, L.; Elboudwarej, O.; Cooke, P. A.; Gray, J. W.; Alivisatos, A. P.; Chen, F. F. Cellular Effect of High Doses of Silica-coated Quantum Dot Profiled with High Throughput Gene Expression Analysis and High Content Cellomics Measurements. Nano Lett. 2006, 6, 800−808. (30) Nagy, A.; Steinbrück, A.; Gao, J.; Doggett, N.; Hollingsworth, J. A.; Iyer, R. Comprehensive Analysis of the Effects of CdSe Quantum Dot Size, Surface Charge, and Functionalization on Primary Human Lung Cells. ACS Nano 2012, 6, 4748−4762. (31) Scialdone, A.; Natarajan, K. N.; Saraiva, L. R.; Proserpio, V.; Teichmann, S. A.; Stegle, O.; Marioni, J. C.; Buettner, F. Computational Assignment of Cell-cycle Stage from Single-cell Transcriptome Data. Methods 2015, 85, 54−61. (32) Darmanis, S.; Sloan, S. A.; Zhang, Y.; Enge, M.; Caneda, C.; Shuer, L. M.; Hayden Gephart, M. G.; Barres, B. A.; Quake, S. R. A Survey of Human Brain Transcriptome Diversity at the Single Cell Level. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, 7285−7290. (33) Achim, K.; Pettit, J. B.; Saraiva, L. R.; Gavriouchkina, D.; Larsson, T.; Arendt, D.; Marioni, J. C. High-throughput Spatial Mapping of Single-cell RNA-seq Data to Tissue of Origin. Nat. Biotechnol. 2015, 33, 503−509.

(34) Zeisel, A.; Munoz-Manchado, A. B.; Codeluppi, S.; Lonnerberg, P.; La Manno, G.; Juréus, A.; Marques, S.; Munguba, H.; He, L.; Betsholtz, C.; Rolny, C.; Castelo-Branco, G.; Hjerling-Leffler, J.; Linnarsson, S. Brain Structure. Cell Types in the Mouse Cortex and Hippocampus Revealed by Single-cell RNA-seq. Science 2015, 347, 1138−1142. (35) Dueck, H.; Khaladkar, M.; Kim, T. K.; Spaethling, J. M.; Francis, C.; Suresh, S.; Fisher, S. A.; Seale, P.; Beck, S. G.; Bartfai, T.; Kuhn, B.; Eberwine, J.; Kim, J. Deep Sequencing Reveals Cell-type-specific Patterns of Single-cell Transcriptome Variation. Genome Biol. 2015, 16, 122−131. (36) Luo, Y.; Coskun, V.; Liang, A.; Yu, J.; Cheng, L.; Ge, W.; Shi, Z.; Zhang, K.; Li, C.; Cui, Y.; Lin, H.; Luo, D.; Wang, J.; Lin, C.; Dai, Z.; Zhu, H.; Zhang, J.; Liu, J.; Liu, H.; deVellis, J.; et al. Single-cell Transcriptome Analyses Reveal Signals to Activate Dormant Neural Stem Cells. Cell 2015, 161, 1175−1186. (37) Kim, D. H.; Marinov, G. K.; Pepke, S.; Singer, Z. S.; He, P.; Williams, B.; Schroth, G. P.; Elowitz, M. B.; Wold, B. J. Single-cell Transcriptome Analysis Reveals Dynamic Changes in lncRNA Expression during Reprogramming. Cell Stem Cell 2015, 16, 88−101. (38) Reyes, J. M.; Chitwood, J. L.; Ross, P. J. RNA-Seq Profiling of Single Bovine Oocyte Transcript Abundance and Its Modulation by Cytoplasmic Polyadenylation. Mol. Reprod. Dev. 2015, 82, 103−114. (39) Piras, V.; Tomita, M.; Selvarajoo, K. Transcriptome-wide Variability in Single Embryonic Development Cells. Sci. Rep. 2014, 4, 7137−7146. (40) Treutlein, B.; Brownfield, D. G.; Wu, A. R.; Neff, N. F.; Mantalas, G. L.; Espinoza, H.; Desai, T. J.; Krasnow, M. A.; Quake, S. R. Reconstructing Lineage Hierarchies of the Distal Lung Epithelium Using Single-cell RNA-seq. Nature 2014, 509, 371−375. (41) Tang, F.; Barbacioru, C.; Bao, S.; Lee, C.; Nordman, E.; Wang, X.; Lao, K.; Surani, M. A. Tracing the Derivation of Embryonic Stem Cells from the Inner Cell Mass by Single-cell RNA-Seq Analysis. Cell Stem Cell 2010, 6, 468−478. (42) Gužvić, M.; Braun, B.; Ganzer, R.; Burger, M.; Nerlich, M.; Winkler, S.; Werner-Klein, M.; Czyż, Z. T.; Polzer, B.; Klein, C. A. Combined Genome and Transcriptome Analysis of Single Disseminated Cancer Cells from Bone Marrow of Prostate Cancer Patients Reveals Unexpected Transcriptomes. Cancer Res. 2014, 74, 7383− 7394. (43) Patel, A. P.; Tirosh, I.; Trombetta, J. J.; Shalek, A. K.; Gillespie, S. M.; Wakimoto, H.; Cahill, D. P.; Nahed, B. V.; Curry, W. T.; Martuza, R. L.; Louis, D. N.; Rozenblatt-Rosen, O.; Suvà, M. L.; Regev, A.; Bernstein, B. E. Single-cell RNA-seq Highlights Intratumoral Heterogeneity in Primary Glioblastoma. Science 2014, 344, 1396−1401. (44) Zhang, X.; Zhang, C.; Li, Z.; Zhong, J.; Weiner, L. P.; Zhong, J. F. Investigating Evolutionary Perspective of Carcinogenesis with Single-cell Transcriptome Analysis. Aizheng 2013, 32, 636−639. (45) Shalek, A. K.; Satija, R.; Shuga, J.; Trombetta, J. J.; Gennert, D.; Lu, D.; Chen, P.; Gertner, R. S.; Gaublomme, J. T.; Yosef, N.; Schwartz, S.; Fowler, B.; Weaver, S.; Wang, J.; Wang, X.; Ding, R.; Raychowdhury, R.; Friedman, N.; Hacohen, N.; Park, H.; May, A. P.; Regev, A. Single-cell RNA-seq Reveals Dynamic Paracrine Control of Cellular Variation. Nature 2014, 510, 363−369. (46) Deng, Q.; Ramskold, D.; Reinius, B.; Sandberg, R. Single-cell RNA-seq Reveals Dynamic, Random Monoallelic Gene Expression in Mammalian Cells. Science 2014, 343, 193−196. (47) Cadima, J.; Jolliffe, I. T. Variable Selection and the Interpretation of Principal Subspaces. Journal of Agricultural, Biological and Environmental Statistics 2001, 6, 62−79. (48) Orr, G.; Panther, D. J.; Phillips, J. L.; Tarasevich, B. J.; Hu, D.; Teeguarden, J. G.; Pounds, J. G. Submicrometer and Nanoscale Inorganic Particles Exploit the Actin Machinery to be Propelled along Microvilli-like Structures into Alveolar Cells. ACS Nano 2007, 1, 463− 475. (49) Xie, Y.; Williams, N. G.; Tolic, A.; Chrisler, W. B.; Teeguarden, J. G.; Maddux, B. L. S.; Pounds, J. G.; Laskin, A.; Orr, G. Aerosolized 10184

DOI: 10.1021/acsnano.6b05452 ACS Nano 2016, 10, 10173−10185

Article

ACS Nano ZnO Nanoparticles Induce Toxicity in Alveolar Type II Epithelial Cells at the Air-Liquid Interface. Toxicol. Sci. 2012, 125, 450−461. (50) Mihai, C.; Chrisler, W.; Xie, Y.; Szymanski, C.; Tolic, A.; Hu, D.; Tarasevich, B.; Orr, G. Intracellular Accumulation Dynamics and Fate of Zinc Ions in Alveolar Epithelial Cells Exposed to Airborne ZnO Nanoparticles at the Air-liquid Interface. Nanotoxicology 2015, 9, 9− 22. (51) Szymanski, C.; Munusamy, P.; Mihai, C.; Xie, Y.; Hu, D.; Gilles, M.; Tyliszczak, T.; Thevuthasan, S.; Baer, D.; Orr, G. Shifts in Oxidation States of Cerium Oxide Nanoparticles Detected Inside Intact Hydrated Cells and Organelles. Biomaterials 2015, 62, 147−154. (52) Ramsköld, D.; Luo, S.; Wang, Y. C.; Li, R.; Deng, Q.; Faridani, O. R.; Daniels, G. A.; Khrebtukova, I.; Loring, J. F.; Laurent, L. C.; Schroth, G. P.; Sandberg, R. Full-length mRNA-Seq from Single-cell Levels of RNA and Individual Circulating Tumor Cells. Nat. Biotechnol. 2012, 30, 777−782. (53) Qiu, S.; Luo, S.; Evgrafov, O.; Li, R.; Schroth, G. P.; Levitt, P.; Knowles, J. A.; Wang, K. Single-neuron RNA-Seq: Technical Feasibility and Reproducibility. Front. Genet. 2012, 3, 124−132. (54) Shalek, A. K.; Satija, R.; Adiconis, X.; Gertner, R. S.; Gaublomme, J. T.; Raychowdhury, R.; Schwartz, S.; Yosef, N.; Malboeuf, C.; Lu, D.; Trombetta, J. J.; Gennert, D.; Gnirke, A.; Goren, A.; Hacohen, N.; Levin, J. Z.; Park, H.; Regev, A. Single-cell Transcriptomics Reveals Bimodality in Expression and Splicing in Immune Cells. Nature 2013, 498, 236−240. (55) Robinson, M. D.; McCarthy, D. J.; Smyth, G. K. edgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data. Bioinformatics 2010, 26, 139−140. (56) Huang, D. W.; Sherman, B. T.; Lempicki, R. A. Systematic and Integrative Analysis of Large Gene Lists using DAVID Bioinformatics Resources. Nat. Protoc. 2009, 4, 44−57.

10185

DOI: 10.1021/acsnano.6b05452 ACS Nano 2016, 10, 10173−10185