Expression Patterns of Energy-Related Genes in Single Cells Uncover Key Isoforms and Enzymes That Gain Priority Under Nanoparticle-Induced Stress.

Cellular responses to nanoparticles (NPs) have been largely studied in cell populations, providing averaged values that often misrepresent the true molecular processes that occur in the individual cell. To understand how a cell redistributes limited molecular resources to achieve optimal response and survival requires single-cell analysis. Here we applied multiplex single molecule-based fluorescence in situ hybridization (fliFISH) to quantify the expression of 10 genes simultaneously in individual intact cells, including glycolysis and glucose transporter genes, which are critical for restoring and maintaining energy balance. We focused on individual gill epithelial cell responses to lithium cobalt oxide (LCO) NPs, which are actively pursued as cathode materials in lithium-ion batteries, raising concerns about their impact on the environment and human health. We found large variabilities in the expression levels of all genes between neighboring cells under the same exposure conditions, from only a few transcripts to over 100 copies in individual cells. Gene expression ratios among the 10 genes in each cell uncovered shifts in favor of genes that play key roles in restoring and maintaining energy balance. Among these genes are isoforms that can secure and increase glycolysis rates more efficiently, as well as genes with multiple cellular functions, in addition to glycolysis, including DNA repair, regulation of gene expression, cell cycle progression, and proliferation. Our study uncovered prioritization of gene expression in individual cells for restoring energy balance under LCO NP exposures. Broadly, our study gained insight into single-cell strategies for redistributing limited resources to achieve optimal response and survival under stress.

Introduction

Lithium cobalt oxide (LCO) is a commonly used cathode material in Li-ion batteries.[1,2] While LCO is used as sintered aggregates of nanoparticles (NPs) with primary particle diameter on the order of 100 nm, under operating conditions these particles fracture into smaller, sheet-like “nanoflakes”.[3] In addition, current efforts are in place to decrease the size of cathode materials to improve battery performance through faster ion and electron transport and increased mechanical stability.[4,5] With little infrastructure in place and low economic incentive to recycle Li-ion batteries,[6] these nanomaterials are likely to end up in landfills, leachate, or as air emissions,[7,8] thus highlighting the need to understand the environmental and health implications of Li-ion battery materials, especially LCO NPs.

It has been shown that exposures of Chironomus riparius, a model organism for aquatic exposures at the benthic zone, to LCO NPs caused a significant decline in larval growth, delay in adult emergence, and reduction in hemoglobin.[9] This study also found significant changes in the expression of genes indicative of metal ion toxicity and oxidative stress, as well as significant reductions in the expression of heme synthesis genes. Further studies in this organism suggested that the oxidation of metabolic and regulatory Fe–S centers within proteins by LCO NPs initiates the disruption of metabolic homeostasis and subsequently the growth and development of the organism.[10] Studies in Daphnia magna, a model organism for freshwater exposures, found that chronic exposures to LCO NPs significantly impacted daphnid reproduction and survival even at low concentrations, as well as dose-dependent downregulations of genes important in metal detoxification, metabolism, and cell maintenance.[11]

Studies in rainbow trout gill epithelial cells, a model cell-type for aquatic environmental exposures, found that LCO NPs led to cell death and strong induction of intracellular reactive oxygen species (ROS).[12] The NPs could be found within abnormal multilamellar bodies and increased the formation of intracellular vacuoles; both are indicative of cellular stress. This study also found a significant increase in p53 gene expression in response to subtoxic doses of the NPs, where no cell death was detected, indicating the presence of stress signals at very low NP concentrations. A followup study aiming to understand the molecular events underlying the strong induction of ROS by LCO NPs in individual cells found a sequence of changes in gene expression, with an initial increase in the expression of genes targeting superoxide species, followed by an increase in the expression of genes targeting peroxide and hydroxyl species, which was consistent with a sequential formation of more complex molecular species.[13]

The health impacts of LCO microscale particles (∼8 μm in diameter) were also studied in mice following lung aspiration,[14] showing acute inflammatory lung responses, which persisted for 2 months in association with fibrosis. Cobalt ions were detected in the broncho-alveolar lavage fluid, associated with upregulation of a marker gene for fibrosis and the biological activity of cobalt ions.[15] Mutagenic activity of these larger particles was found in mice and lung cells in cultures,[16] where hydroxyl radicals, DNA strand breaks, and oxidative lesions were induced, likely by the cobalt ions that were released from these larger particles. Interestingly, studies of nanoscale LCO particles in the environmental model systems mentioned above, including Chironomus,[9]Daphnia,[11] and gill epithelial cells,[12,13] showed that exposures to lithium or cobalt ions alone had no impact on the viability of the organisms or cells, as well as on ROS generation, pointing to the intact LCO NPs as the cause for these adverse effects.

Recently, we conducted cell population transcriptomic analyses (bulk RNA-Seq) of trout gill epithelial cells exposed to LCO NPs at a subtoxic dose (showing no impact on cell death) and a toxic dose (the lowest dose to induce cell death) for 24 and 48 h, as well as to lithium and cobalt ions at concentrations released from the toxic dose.[17] Only a few genes were found to be impacted by the ions at 24 h but were restored to normal levels at 48 h. In contrast, strong upregulation of energy-related processes, such as glycolysis and glucose transport processes, as well oxygen and hypoxia-related processes, was observed in response to both toxic and subtoxic NP doses. These observations point to cellular attempts to restore oxygen and energy imbalance by increasing glycolysis rate and glucose import and transport.[17,18] However, how the cells prioritize the expression of these genes to achieve optimal cellular response and survival with limited resources under the NP insult is still unclear. More broadly, how a cell under critical stress conditions redistributes resources to restore and maintain energy balance is unknown.

While cell population analyses, such as bulk RNA-Seq, could provide averaged values and general trends, such approaches cannot uncover accurate comparisons between the expression levels of the different genes or prioritization of molecular resources as they truly occur in individual cells. This limitation of bulk analyses is mainly due to the large variability between cells, especially when different cell subpopulations are present. For example, if one subpopulation upregulates the expression of gene A and downregulates the expression of gene B, while the other subpopulation upregulates the expression of gene B and downregulates the expression of gene A, the averaged values from bulk RNA-Seq can show little or no change in the expression of both genes. The patterns of gene expression in individual cells usually span a wide range of possibilities leading to complex relationships or ratios between the expression of the genes.

The presence of multiple cell subpopulations is especially relevant to cells exposed to NPs. It has been shown that NP exposures end with a wide range of NP loads in individual cells, where some cells might carry no or only a few NPs, while other cells can be loaded with hundreds or even thousands of NPs per cell.[19−21] Such wide distribution in NP load per cell has been shown to correlate with a wide distribution in the level of responses.

To understand how the cell optimizes gene expression to restore and maintain energy balance with limited resources under the NP insult, we quantified simultaneously the expression of 8 glycolysis and 2 glucose transporter genes in intact gill epithelial cells using fluctuation localization imaging-based fluorescence in situ hybridization (fliFISH), a highly accurate single molecule-based FISH approach.[13,22] While this approach is limited in the number of genes that it can quantify simultaneously in individual cells, its advantage, compared to sequencing-based approaches, is the ability to accurately detect very low abundance genes, down to single transcripts, with no cell destruction or amplification biases.

We found large variabilities in the expression levels of all 10 genes between neighboring cells in the dish under all exposure conditions, which was manifested in single cell clustering analysis, where exposed cell subpopulations were found both among and separate from untreated control cells. By quantifying gene expression ratios between pairs of genes in individual cells, we identified shifts in favor of genes with key roles in the cellular attempt to restore and maintain energy balance under the NP insult. These prioritized genes include isoforms that can secure or increase glycolysis rate more efficiently, as well as genes with multiple cellular functions in addition to their function in glycolysis, such as DNA repair, regulation of gene expression, and cell cycle progression and proliferation, among other processes. Thus, the shifts in expression ratios highlighted genes that play critical roles under NP-induced stress, pointing to key isoforms and enzymes with multiple cellular roles that gained priority for cell survival under stress.

Results and Discussion

Nanoparticle Characterization

Here we used LCO NPs from the same NP batch used in Mensch et al.,[17] where a detailed characterization of the NPs is provided. Briefly, the NPs used here have stoichiometry of Li0.87CoO2, as determined by ICP-OES. XPS spectra for the Co 2p and Li 1s regions confirmed the chemical composition of the LCO NPs. These NPs have a sheet-like morphology as determined by TEM and SEM, and a thickness of ∼5 nm in both Nanopure H2O and growth medium, as determined by AFM. The ζ-potential of the particles is −7 ± 1 mV in Nanopure water and −10 ± 1 mV in growth medium.

Large Variabilities in the Expression Levels of All 10 Genes Were Found between Individual Cells under the Same Exposure Conditions

In this study, we used a high-accuracy single-molecule-based fluorescence in situ hybridization (FISH) approach, fluctuation localization imaging-based FISH (fliFISH),[13,22,25] to quantify the expression of selected genes in intact rainbow trout gill epithelial cells—a model cell-type for aquatic environmental exposures.[23,24] As described in Figure , we used 5 fluorescent colors to create 10 distinct two-color barcodes to target each of the 10 selected genes in each cell. Eight of the genes encode enzymes in the glycolysis pathway and 2 genes encode glucose transporters, together participating in key processes for producing and maintaining energy balance. The 10 genes were selected based on previous observations using bulk RNA-seq in these cells, showing a substantial increase in their expression levels in response to LCO NPs.[17] The full names for the abbreviated genes and the key for the color-codes assigned to each gene through the manuscript are provided in Figure .

Figure 1

Combinatorial barcoding fliFISH is used to count transcripts of glycolysis and glucose transporter genes in individual cells. A. Five fluorescent dyes (colored circles) are used to create 10 distinct two-color barcodes to target the selected 10 genes (colored boxes). Each gene is then assigned with a distinct pseudocolor (color-code), which is used throughout the manuscript. B. An example for a fliFISH image where individual cells are outlined with the dashed line and the nuclei are stained in blue. Transcript molecules for each gene are detected by the colocalization of their respective two colors (dots) with nanometer resolution. The number of transcripts are then counted in individual cells providing the expression level for each gene. The area marked by the square in the left image is enlarged in the right image. Color-Code Key: glucose transporter 1 (GLUT1); glucose transporter 4 (GLUT4); hexokinase-1 (HK1); hexokinase-2 (HK2); phosphofructokinase 1 (PFK1); phosphofructokinase 2 (PFK2); fructose-bisphosphate aldolase A (ALDOA); glyceraldehyde-3-phosphate dehydrogenase (GAPDH); phosphoglycerate mutase 1 (PGAM1); pyruvate kinase (PKM).

Here we studied the responses of the cells to LCO NPs at a subtoxic dose (5 μg/mL), showing no impact on cell viability, and a toxic dose (25 μg/mL), which was the lowest dose to induce significant cell death, at both 24 and 48 h (Figure S1). An example for the distribution of single cell gene expression levels, determined using fliFISH, is provided in Figure and Figure S2. Each cell is assigned a color intensity or shade in the respective gene color based on the number of transcripts found in the cell, normalized to the highest single cell transcript counts found under all conditions for that gene. Thus, cells with relatively high transcript counts are shown in intense or bright colors, while cells with relatively low transcript counts are shown in faint or light colors. These maps demonstrate the large variability that we found in the expression levels of all 10 genes between neighboring cells under each exposure condition, from only a few copies to over a hundred copies in individual cells for certain genes (Figures and S2). Interestingly, these maps show that individual cells often express high levels of all 10 genes, while other cells express low levels of the 10 genes. Rarely does a cell express high levels of some genes and low levels of other genes. This consistency in single cell expression levels across the 10 genes is observed under all exposure conditions, including in control, unexposed cells.

Figure 2

Mapping gene expression levels in single cells. Each series, consisting of 10 images, shows the same cells, colored with the respective gene color-code (key on the left). Each cell was assigned an intensity level or shade of the respective gene color based on the number of transcripts found in the cell (number placed on each cell), normalized to the highest single cell transcript counts found under all conditions for that gene. Thus, cells with high transcript counts are shown in bright colors, while cells with low transcript counts are shown in light colors. Scale bars equal 10 μm. The upper, middle, and lower series show examples for control cells, cells exposed to the subtoxic dose for 24 h, and cells exposed to the toxic dose for 24 h, respectively. Series for cells exposed to the subtoxic and toxic doses for 48 h are shown in Figure S1.

However, when averaged over ∼100 cells per exposure condition, the averaged single cell gene expression levels showed strong upregulation of all 10 genes in cell exposed to the toxic dose (25 μg/mL) for 24 and 48 h compared to control cells (Figure S3). A significant increase in the averaged expression level of a subset of the genes was also observed in response to the subtoxic dose (5 μg/mL) at 48 h, with minimal changes detected at 24 h, together pointing to an evolution in the response to this subtoxic dose over time (Figure S3).

Consistent with the variability observed in the expression levels between neighboring cells (Figures and S2), clustering analysis of gene expression in individual cells showed a wide distribution under all exposure conditions, where some cells clustered together with control cells, while other cells clustered away from the control cells (Figure A–D). In Figure plots, individual dots represent individual cells (blue dots for control cells and red dots for exposed cells), where the distances between dots or cells are based on cell size-normalized gene expression of the 10 genes.

Figure 3

Clustering analysis by single cell raw gene expression levels. Each cell is represented by a dot—blue dots for control cells and red dots for exposed cells. Distance between individual cells (dots) is based on cell size-normalized gene expression of the 10 genes. Distance is calculated using the Euclidian calculation in a multidimensional scaling plot. Ellipses are designed to encompass 95% of all dots for each condition.

The ellipses encompass 95% of control (blue ellipses) and exposed (red ellipses) cells. This analysis clearly shows that some of the exposed cells had distinct response patterns while other cells had similar response patterns as control cells under all exposure conditions. One way to explain such distribution in the exposed cells is based on the previously reported heterogeneity of the NP load in individual cells, where some cells might carry no or only few NPs, while other cells can be loaded with hundreds or even thousands of NPs.[19−21] Such a wide distribution in NP load per cell has been shown to correlate with a wide distribution in the level of responses.

To better understand the distinct responses of the exposed cells compared to control cells, we narrowed the control ellipses to encompass 40% of the control cells while keeping them at the center of the distributions to represent the core response of the control cells (Figure A–D). Comparing responses of exposed cells outside the ellipses to control cells inside the ellipses could shed light on how distinct subpopulations of individual cells prioritize gene expression to optimize energy balance and survival under the NP insult. The averaged change in gene expression levels in exposed cells outside the ellipses compared to control cells inside the ellipses (Figure ) showed similar patterns to those generated by comparing all exposed to all control cells (Figure S3), but resulted in noticeably more robust changes in responses. Interestingly, the exposed cells that overlapped with the control cells inside the ellipses also showed changes in the averaged expression levels of certain genes compared to the control cells, mainly downregulation (Figure S4), although not nearly as significant as the changes observed in the exposed cells outside the ellipses.

Figure 4

Clustering analysis where distance between individual cells (dots) is based on cell size-normalized raw gene expression of the 10 genes. Distance is calculated using the Euclidian calculation in a multidimensional scaling plot. Blue dots represent control cells and red dots represent treated cells. The blue ellipses are designed to encompass 40% of the control cells while being kept at the center of the distributions to represent the core response of the control cells. Comparing responses of exposed cells outside the ellipse to responses of control cells inside the ellipse will help understand how distinct subpopulations of cells prioritize gene expression to restore energy balance and survive under the NP insult.

Figure 5

Fold change (A) and significance (B) in gene expression levels for the 10 genes (color-coded), averaged across all the exposed cells outside the ellipse, compared to the control cells inside the ellipse. The exposure conditions include a toxic dose (25 μg/mL) for 24 and 48 h, and a subtoxic dose (5 μg/mL) for 24 and 48 h. Significance is represented as the −log10 of the p-value of the difference, and the horizontal line indicates significance (p-value = 0.05).

Shifts in Expression Ratios between Gene-Pairs in Individual Cells Identified Enzymes That Gained Priority under NP-Induced Stress

To understand how cellular resources are allocated within an individual cell to keep glycolysis and glucose transport functioning under NP-induced stress, we first plotted the raw expression levels of the 10 genes in individual cells for each of the 10 genes (Figure (A,B). Each bar in the graphs of Figure represents a cell, where bars on the left of the black line represent control cells inside the ellipses, while bars on the right represent exposed cells outside the ellipses. The patterns that emerge point to shifts in gene expression levels and resulting ratios that occur in the exposed cells compared to control cells. To quantify such shifts, we calculated raw gene expression ratios of all possible pairs between the 10 genes in each cell (total of 45 pairs per cell) and assessed the degree of changes in these ratios in exposed cells (outside the ellipses) compared to control cells (inside the ellipses). The top 4 most significantly changed gene pair ratios for each exposure group compared to the control are shown in Figure A–D, where each gene expression level in the pair is displayed as a fraction of their total level. In other words, we zoomed-in on gene pairs in Figure whose ratios changed the most in the exposed cells compared to control cells. This analysis allowed us to understand how the cell prioritizes the expression of glycolysis and glucose transport genes for optimal response and survival under the insult. While multiple highly changed ratios could be detected as a result of the exceptionally high or low expression of a single gene, focusing on the 4 most significantly changed ratios allowed us to identify and focus on relationships between genes that best highlighted the cellular priorities among these 10 genes.

Figure 6

Raw expression levels of the 10 genes in individual cells, where each bar in the graphs represents a cell. Bars on the left of the black line represent control cells inside the ellipse, while bars on the right represent exposed cells outside the ellipse. A. Gene expression levels in response to the toxic dose at 24 h (left) and 48 h (right). B. Gene expression levels in response to the subtoxic dose at 24 h (left) and 48 h (right).

Figure 7

Comparing changes in gene expression ratios within individual cells. Raw gene expression ratios were calculated between all possible pairs among the 10 genes (45 pairs) in each cell. Differences in the ratios for each of these gene pairs was then assessed in exposed cells (outside the ellipses shown in Figure ) compared to control cells (inside the ellipses in Figure ). The top 4 most significantly changed gene pair ratios for each exposure condition compared to control are shown. Each bar represents a cell. The expression level of each gene in the pair is displayed as a fraction of their total level in the cell. Bars representing control cells (inside the ellipse) are shown on the left of the vertical black line, while bars representing exposed cells (outside the ellipse) are shown on the right. The averaged expression levels for each pair are plotted on the left of each graph and marked as “C” for control cells and “E” for exposed cells.

HK1 and HK2

One of the most significantly changed ratios was observed between the expression levels of hexokinase-2 (HK2) and hexokinase-1 (HK1). Figure A,B,D shows a significant shift in favor of HK1 expression (green) relative to HK2 expression (magenta) in response to all exposure conditions except the subtoxic dose at 24, suggesting an increased priority of HK1 relative to HK2 expressions under NP induced stress. This interpretation is also supported by the averaged expression levels of the two genes presented in Figure , showing a greater increase in HK1 under all exposure conditions.

Biological Meaning in the Context of the Literature

Hexokinases (HKs) are a family of 4 isoforms that, among other functions, catalyze the first step in glycolysis by converting glucose to glucose 6-phosphate (G-6-P). The two most common isoforms are HK1 and HK2, which in mammalian systems are expressed in overlapping tissues but have different subcellular locations.[26,27] HK1 is mainly associated with mitochondria, while HK2 is associated with mitochondria and other cytoplasmic compartments. It has been shown that the different subcellular distributions are associated with different metabolic roles, where mitochondrial-bound HKs channel G-6-P toward glycolysis (catabolic use), while cytoplasmic HKs regulate glycogen formation (anabolic use).[28] It has been suggested that unlike HK2, which contributes to both catabolic and anabolic functions, the mitochondrial-bound HK1 is mainly committed to glycolysis, together enabling cells to adapt to changing metabolic conditions while maintaining energy balance.[26] It is possible that the increase in the expression of HK1 relative to HK2 expression in cells exposed to the NPs reflects an attempt to secure and restore energy balance by prioritizing glycolysis or direct energy production over other cellular processes.

PFK1 and PFK2

Among the 4 most significantly shifted ratios was the ratio between the expressions of phosphofructokinase 1 (PFK1) and phosphofructokinase 2 (PFK2). Figure A,B shows that in control cells, on average, the expression of PFK1 (turquoise) was higher compared to PFK2 (periwinkle). However, this ratio was significantly shifted in favor of PFK2 in cells exposed to the toxic dose at both time points, bringing the expressions of both PFK2 and PFK1 to nearly equal levels at 48 h. Although the averaged expression level of both PFK1 and PFK2 increased in response to both toxic and subtoxic doses (Figure ), a significant shift in favor of PFK2 relative to PFK1 was detected mainly in response to the toxic dose, suggesting an increased priority of PFK2 expression under a strong NP induced stress.

PFK1 is a key regulator of the overall glycolysis reactions. It catalyzes the first distinct step in glycolysis, converting fructose-6-phosphate to fructose-1,6-bisphosphate. The second isoform, PFK2, catalyzes the conversion of fructose-6-phosphate to fructose-2,6-bisphosphate, which is a stimulator/activator of PFK1.[29,30] When PFK2 is phosphorylated, it acts as a phosphatase and decreases the concentration of fructose-2,6-bisphosphate by converting it back to fructose-6-phosphate, resulting in inhibition of glycolysis and stimulation of gluconeogenesis. When PFK2 is nonphosphorylated, it acts as a kinase and increases the level of fructose-2,6-bisphosphate, which in turn increases PFK1 activation and subsequently stimulation of glycolysis.[31,32] It is possible that the significant shift in favor of PFK2 expression relative to PFK1, detected in cells exposed to the toxic dose, reflects an attempt to stimulate glycolysis via an increase in fructose-2,6-bisphosphate generation and activation of PFK1. The PFK2 family of enzymes consists of 4 members (PFKFB1–4), with PFKFB3 being expressed in our system. In addition to playing an important role in glycolysis, PFKFB3 has been shown to play a key role in cell cycle regulation and prevention of apoptosis,[33,34] as well as in supporting cell proliferation and migration in several systems, including zebrafish embryos.[35] It is possible that the significant shift in favor of PFK2 (PFKB3) expression relative to PFK1 in cells exposed to the toxic dose at both time point (Figure A,B), reflects the additional role of PFK2 in regulating cell cycle and supporting cell proliferation under the NP insult.

GAPDH

A significant shift was also observed in favor of glyceraldehyde-3-phosphate dehydrogenase (GAPDH) (cyan) in cells exposed to the toxic dose at both time points compared to control cells (Figure A,B).

Biological Meaning in the Context of the Literature

GAPDH is often described as a “jack-of-all trades”.[36] In addition to playing an important role in glycolysis, GAPDH localizes to multiple cellular compartments and plays key roles in many other cellular processes.[37] GAPDH plays key roles in DNA replication and repair,[38] regulation of gene expression,[39] cell cycle progression,[40] redox sensing,[41] and apoptosis,[42] among other processes. These functions are likely activated as part of the cellular response to the NP insult, which has been reported to include strong ROS generation and oxidative stress, DNA and heme protein damage, and energy and oxygen related imbalance.[10,12,16,17] This interpretation is also supported by the observation that the averaged expression level of GAPDH significantly increased in response to all exposure conditions except the nontoxic dose at 24 h (Figure ).

PGAM1

Among the topmost significantly changed ratios was a shift in favor of phosphoglycerate mutase 1 (PGAM1) (olive green) when paired with three other genes, including HK2 (Figure B,C,D), fructose-bisphosphate aldolase A (ALDOA), and glucose transporter 4 (Glut4) (Figure B,D). On average, PGAM1 expression levels significantly increased in cells under all exposure conditions except under the subtoxic dose at 24 h (Figure ).

Biological Meaning in the Context of the Literature

PGAM1 has been found to be overexpressed in multiple cancer tissues and promote cell proliferation and cancer progression, thus being a potential therapeutic target.[43,44] It has been found that both the glycolytic and nonglycolytic functions of PGAM1 underlie its role in promoting cell proliferation and cancer progression.[45] Our observed increase in the priority of PGAM1 expression in cells under NP-induced stress might reflect the key role of PGAM1 in increasing the rate of glycolysis required for cell proliferation.

GLUT1 and GLUT4

The averaged expression levels of both glucose transporters 1 (GLUT1) and 4 (GLUT4) increased significantly in response to the toxic dose at both time points, but only GLUT1 showed significant increase in response to the subtoxic dose at 48 h (Figure ). Interestingly, two of the four most significantly changed ratios in cells exposed to the subtoxic dose at 48 h showed decrease in GLUT4 (orange) expression relative to HK1 and PGAM1 (Figure D), which have been implicated in securing or increasing glycolysis rate, as described above.

Biological Meaning in the Context of the Literature

Both GLUT1 and GLUT4 facilitate the passive movement of glucose down the concentration gradient.[46] GLUT4 is present mainly in insulin-sensitive tissues, where insulin stimulation increases GLUT4 translocation from intracellular compartments to the cell membrane.[47] In contrast, GLUT1 is present in most cells, where it is responsible for the basal glucose uptake.[48] It has been shown in mammalian systems that several stress stimuli, such as hypoxia[49] and osmotic shock[50] increase GLUT1 levels. It has been also shown that under certain conditions, GLUT1 and GLUT4 expression is regulated in opposite directions.[51] Prior observations suggesting that LCO NPs induce hypoxia and other stress responses in these cells[17] might explain the more consistent increase in the expression of GLUT1 (Figure ) while decreasing the priority of GLUT4 relative to other more critical genes for cell survival under the NP insult (Figure D).

The expression level of actin was quantified in a different set of single cells exposed to toxic and subtoxic NP doses for 24 and 48 h, as well as in control, unexposed cells (Figure S5). As expected, the results show, on average, comparably high expression levels of actin across all conditions. To gain additional insights to the physiological state of these cells, the expression of MKI67, encoding the proliferation marker, Ki-67, was also quantified in the same single cells (Figure S5). In contrast to the high levels of actin, MKI67 showed lower expression levels, and on average, a significant downregulation in response to both the toxic and subtoxic doses at 24 h. Such downregulation indicates a shift of resources away from proliferation while giving priority to genes supporting cellular recovery. At 48 h, MKI67 showed normal levels, potentially reflecting some regaining of cellular functions, at least in a subset of cells that might have been less loaded or impacted by the nanoparticles. In an earlier study focusing on the impact of the same LCO NPs in the same cell line, it was found that the expression of P53, involved in regulating DNA repair and apoptosis, was upregulated even in response to a subtoxic dose of 1 μg/mL.[12] Together, these results indicate shifting resources away from proliferation while increasing priority of genes responsible for repairing and restoring cellular functions.

Conclusion

One of the more important findings in our study is the identification of key genes whose expression gained priority over other genes in the same cell under LCO NP-induced stress. This information, which requires single-cell analysis, sheds light on the role of these genes in restoring energy balance under LCO NP exposure and, more broadly, uncovers single-cell strategies for redistributing limited resources to achieve optimal response and survival under stress. By quantitative single-cell gene expression analysis, we found that the genes that gain priority in response to LCO NP-induced stress are either more efficient in contributing to restoring energy balance or play important roles in multiple cellular functions that are critical for cell survival.

We found large variabilities in the expression levels of all 10 genes between neighboring cells under all exposure conditions, from only a few copies to over 100 copies in individual cells for certain genes. Interestingly, individual cells often expressed all 10 genes at high levels or all 10 genes at low levels, rarely showing high levels of some genes and low levels of other genes in the same cell. Consistent with the variability observed between neighboring cells in the dish, clustering analysis of gene expression in individual cells showed wide distributions under all exposure conditions, with some cells clustering together with control cells, while other cells showed distinct expression patterns. Such variability could be explained by the previously reported heterogeneity of NP load in individual cells.

Averaging single cell expression values over ∼100 cells per exposure condition showed robust upregulations of all 10 genes in response to the toxic dose at both time points, pointing to strong energy imbalance. Significant increases in the averaged expression level of a gene subset was also observed in response to the subtoxic dose at 48 h, while minimal change was detected at 24 h, pointing to an evolution in the response to this subtoxic dose over time.

Comparing gene expression ratios between all possible gene pairs among the 10 genes in each cell uncovered shifts in favor of genes that play key roles in the cellular attempt to restore and maintain energy balance under the NP insults. Among these genes are HK1 and PFK2, which gained priority relative to their isoforms, HK2 and PFK1, respectively, likely due to their greater roles in securing and increasing the rate of glycolysis. PFK2 potentially contributed also via its additional roles in regulating cell cycle and supporting cell proliferation. Two other genes gained priority likely due to their other key roles, in addition to their role in glycolysis. These include PGAM1, likely reflecting its roles in promoting cell proliferation, and GAPDH, likely reflecting its roles in DNA repair, regulation of gene expression, and cell cycle progression, among other processes. Thus, the shifts in expression ratios highlighted genes that play critical roles under NP-induced stress, pointing to more specialized isoforms or enzymes with multiple cellular roles that gained priority for cell survival under stress.

Methods

NP Characterization

In the current study, the same batch of LCO NPs was used as in Mensch et al.[17] where a detailed characterization of the NPs is provided using SEM, TEM, XPS, AFM, and LDM. These NPs were synthesized as described previously.[52]

LCO NP Preparation for Cell Exposures

A 1 mg·mL–1 stock solution was prepared in growth medium. NP solutions were sonicated in ice water using a Misonix Sonicator 3000, operated at 10 W for 4 × 2.5 min. Stock solution was immediately diluted into sonicated growth medium for cell exposures.

Cell Culture Growth and Exposure

Oncorhynchus mykiss (rainbow trout) gill epithelial cells (RTgill-W1, ATCC CRL-2523) were cultured in Leibovitz’s L-15 growth medium (ATCC) supplemented with 1% antibiotics and 10% fetal bovine serum (ThermoFisher), referred to as “growth medium” through the manuscript. Cells were incubated in ambient atmosphere at 19 °C. Cells were seeded in 35 mm glass coverslip bottom dishes (P35G-1.5-20-C, MatTek life science) until they reached near 90% confluency, when they were exposed to the sonicated NPs suspensions at 5 μg/mL or 25 μg/mL for 24 and 48 h. Cells where then fixed with 4% paraformaldehyde (PFA) and further processed for FISH hybridization. Four plates were used for each exposure condition, including for control unexposed cells.

fliFISH Concept

fliFISH takes advantage of photoswitchable dyes and super-resolution localization microscopy to accurately count and localize mRNA molecules using a small number of oligonucleotide probes.[13,22] The single-molecule on-time fraction (duty cycle) of the fluorescent dyes is measured using optimized excitation conditions.[13] Following probe hybridization and imaging, distinct photoblinking patterns or ensemble on-time fractions are detected from fluorescent spots. These ensemble on-time fractions can distinguish true signals from background noise. True signals, coming from hybridized probes, will show the expected ensemble blinking patterns, estimated from the average on-time fraction of a single probe multiplied by the number of probes used to target a transcript. In contrast, noise from stray or nonspecifically bound probes would generate near single-molecule on-time fraction values or less, and autofluorescence or aggregated probes would rarely generate blinking patterns. Because fliFISH is imaged using the STORM technique, it can also resolve multiple transcripts in a diffraction-limited area.

fliFISH Probe Design

fliFISH was used following the approach describe earlier.[13,22] Each primary FISH probe contained a sequence of ∼20 oligonucleotides (NTs) complementary to the target mRNAs (target domains), extended on each side by two different sequences of 28 NT overhangs complementary to the secondary probes. Eight primary probes were designed to target each gene’s mRNA. The target domain sequences of all probes used in this study are provided in Table S1. Secondary probes, each tagged with two photoswitching dye molecules of the same color (Atto-488, TMR, Alex-594, Alex-647, or Alexa-750), were used to hybridize with the 5′ or 3′ overhang sequences to generate 2 color barcodes (see illustration in Figure ). The sequences of the secondary probes are provided in Table S2. All probe sequences were subjected to BLAST searching to avoid nonspecific targeting and purchased from Integrated DNA Technologies. The hybridization protocol has been described in detail previously.[13,25,53] Briefly, primary probes were first hybridized with secondary probes in a tube to form fluorescent complexes and then introduced to the cells. Following fixation in ice cold 4% PFA, the cells were permeabilized with 70% ethanol. 60 nM probe complexes were used to hybridize with the cells at 37 °C overnight in a humid chamber. The hybridization solution contained 1× SSC, 15% formamide, 10% dextran sulfate, 2 mM ribonucleoside vanadyl complex, 3.4 mg/mL tRNA, 0.2 mg/mL RNase-free BSA. The next day, cells were thoroughly washed with 15% formamide in 1× SSC, and counter-stained with DAPI.

Fluorescence Microscopy

Imaging was performed using Olympus IX-71 inverted microscope with a 100× oil immersion objective, six solid-state lasers (405, 488, 542, 594, 640, 730 nm), and an Electron Multiplying CCD camera (Andor iXon Ultra 897). The camera’s pixel size was 16 μm, and the microscope magnification was 100×. Thus, the scale factor was 160 nm/pixel. 2000 images of 512 × 512 pixels were collected within 400 s and stored as 16-bit FITS files. Imaging was done in an oxygen-depleting buffer described earlier.[22] Aberration between color channels was corrected by calibration of the system with broadband emission fluorescence beads.

fliFISH Data Analysis

The method for identifying and quantifying transcripts in individual cells has been described in detail previously.[13,25,53] Briefly, the centroid of each photoswitching fluorescent spot showing on–off emission or a blinking pattern was determined from the imaging stack in each of the 5 fluorescence channels corresponding to the 5 fluorescent dyes. This was done using a Gaussian musk fitting algorithm to find the central location of each emission event. Blinking density maps of nearby events were then generated, where the density value of each pixel indicates the number of blinking events within a distance R of that pixel. Such a map was used to group blinking events into clusters. R was determined by the length of the target mRNAs and the probe localization error upon hybridization, which in our study was determined to be ∼30 nm. The center of mass in each cluster was evaluated for potential representation of a transcript or multiple transcripts. To confirm the presence of mRNA molecules, two centers of mass from the corresponding two-color channels for each gene had to colocalize within an area of 2× localization errors (∼60 nm). The analysis was performed using MATLAB routines available upon request.

Bioinformatics

Single cell transcript counts, acquired by fliFISH, were normalized according to cell size by dividing the counts for each of the 10 genes by the measured area of each cell’s imaged cross section. This essentially converted transcript counts in each cell to an estimation of single cell transcript concentration. Tests for statistical difference between single cell subpopulations were assessed using the t.test() base function in the R statistical programming language (referred to as R).[54] Bar graphs and scatterplots were generated using the ggplot2 package for R.[55]

Placement for individual cells was determined using the base R function dist(), which by default uses a Euclidian distance approach, and then applying a multidimensional scaling (MDS) approach using the base R function cmdscale(). MDS plots use Euclidian distance as input and assign each sample to a point on a Cartesian grid so that the distances between samples are preserved as much as possible. Ellipses were added using the ggplot2 function stat_ellipse() with parameters type = ”norm” and level = 0.95 (Figure ) or 0.4 (Figure ). This function draws an ellipse over the data in such a way that a specified fraction of the central-most points in the data set are encompassed by the ellipse, in our case, 95% and 40% in Figures and 5, respectively. We used this approach to identify the “core” state of control or treated cells. Cells were identified as inside or outside of ellipses using the function point.in.polygon() from the R package sp.[56,57]

To determine the top 4 gene ratio changes for each condition, the ratio for each gene pair in all treated cells was compared to comparable ratios for all control cells. The top 4 most significantly changed ratios, as determined using the t.test() function above, were selected for further evaluation. For example, if the ratios of gene 1 to gene 2 in control cells are on average close to 1 (i.e., the gene expression levels are similar), but the ratios of gene 1 to gene 2 in treated cells are on average close to 3, this would indicate a noticeable change in the ratio of gene 1 to gene 2 with treatment, and this ratio would be likely to be flagged as a top gene ratio change in our assessment. For Figures and 7, gene levels are depicted as the fraction of the total level of all 10 genes (Figure ) or the current pair (Figure ). The black lines showing the separation between control and treated cells do not line up in Figure mainly because different numbers of exposed cells are compressed into the same figure width, and also because the ellipses, designed to encompass 40% of control cells, do not always include a completely identical set due to differences in cell layout in each plot of the four exposure conditions.