An Integrated Preprocessing Approach for Exploring Single-Cell Gene Expression in Rare Cells

Exploring the variability in gene expressions of rare cells at the single-cell level is critical for understanding mechanisms of differentiation in tissue function and development as well as for disease diagnostics and cancer treatment. Such studies, however, have been hindered by major difficulties in tracking the identity of individual cells. We present an approach that combines single-cell picking, lysing, reverse transcription and digital polymerase chain reaction to enable the isolation, tracking and gene expression analysis of rare cells. The approach utilizes a photocleavage bead-based microfluidic device to synthesize and deliver stable cDNA for downstream gene expression analysis, thereby allowing chip-based integration of multiple reactions and facilitating the minimization of sample loss or contamination. The utility of the approach was demonstrated with QuantStudio digital PCR by analyzing the radiation and bystander effect on individual IMR90 human lung fibroblasts. Expression levels of the Cyclin-dependent kinase inhibitor 1a (CDKN1A), Growth/differentiation factor 15 (GDF15), and Prostaglandin-endoperoxide synthase 2 (PTGS2) genes, previously shown to have different responses to direct and bystander irradiation, were measured across individual control, microbeam-irradiated or bystander IMR90 cells. In addition to the confirmation of accurate tracking of cell treatments through the system and efficient analysis of single-cell responses, the results enable comparison of activation levels of different genes and provide insight into signaling pathways within individual cells.


Results
We first performed on-chip RNA capture validation, then characterized on-chip RT efficiency and bead-based reaction efficiency, as well as quantifying the releasing efficiency of cDNA from beads. We then demonstrated single-cell analysis with the integrated system from microbeam irradiation, cell picking, and on-chip reactions, all the way through QuantStudio dPCR. The integration of the system allows the exploration of cell-to-cell variability in response to microbeam irradiation and in bystander cells.
Validation of RNA capture on the Pre-GEM microchip. Before characterizing efficiencies of downstream reactions, the on-chip capture of mRNA via 5′-PC-Biotin-(dT) 25 -3′ Streptavidin magnetic beads was validated using XenoRNA template. The value of ΔR n , indicating the magnitude of the fluorescent signal and, therefore, amplification generated by PCR, was 3.7 for the positive control (10 5 XenoRNA with 3.5 × 10 6 beads) after 40 cycles of bead-based PCR. Referring to the method section, the remaining solutions after bead capture were mixed with bead-based solution, and RT-PCRs were performed to quantify remaining templates in the On-chip RT characterization and interference of beads with reaction efficiency. Starting with 20k, 50k and 100k copies of XenoRNA, the on-chip bead-based RTs followed by off-chip bead-based qPCRs were conducted to test repeatability and characterize reaction efficiency (Fig. 2a). We found that under the given experimental conditions, the PCR efficiency defined by (10 −1/k −1) × 100%, where k is the slope of the C t as a function of the logarithm of the template copy number, for our bead-based approach was 90.5%. To characterize the on-chip RT efficiency, we conducted the bead-based PCR testing following in-tube RT. The PCR efficiency E b for the bead-based PCR with in-tube RT was 83.2%, meaning the bead-based on-chip RT has a 1.08-fold increase in efficiency compared with that of in-tube RT. This improved efficiency for on-chip reaction likely resulted from more efficient molecular interactions in the microscale environment 25,40 . Meanwhile, we ran a solution-based PCR without introducing beads following in-tube RT. Some bead interference was observed, affecting the reaction efficiency, as the presence of beads during the PCR reaction yielded higher Ct values than the solution-based reactions. In terms of the reaction efficiency, it was found in Fig. 2b that under the given experimental conditions, the PCR efficiency E s for the solution-based PCR testing (86.3%) was higher than the PCR efficiency E b for the bead-based PCR (83.2%), which was likely attributable to enhanced molecular interactions in the absence of the beads. The reaction efficiency was not dependent on template copy number, however. Based on these results, we calculated a correction factor that can be used to estimate the corrected threshold cycle of the bead-based approach in order to offset the influence of beads on the reaction. The corrected threshold cycle can be evaluated by Ct s = Ct b × log(1 + E b )/log(1 + E s ), where Ct b , Ct s , E s , and E b denote threshold cycle of the bead-based positive controls, threshold cycle of the corrected positive controls, solution-based PCR efficiency, and bead-based PCR efficiency, respectively.  www.nature.com/scientificreports www.nature.com/scientificreports/ Characterization of the releasing efficiency of cDNA from the beads via on-chip photocleavage. It is also important to know the releasing efficiency of cDNA from the beads via on-chip photocleavage.
Therefore, an experiment was designed to characterize the release efficiency (Fig. 3). The first arm of the experiment used on-chip cell lysis, mRNA capture, and RT processing of a small number (10) of IMR-90 cells obtained by dilution. Without photocleavage, the RT products (cDNA on beads) were collected and analyzed by qPCR. We next performed on-chip lysis of 10 cells, mRNA capture and RT, followed by on-chip photocleavage via UV irradiation (peak wavelength: 365 nm). The beads were separated from the cleaved RT product (cDNA), and qPCR analyses were performed separately using the cDNA cleaved from the beads and the remaining beads. Looking at the amplification curves for the housekeeping gene GAPDH, there is negligible amplification of the cleaved beads, indicating the cDNA has been effectively removed from the beads by the UV irradiation. We have also noticed that the threshold cycle Ct b of cDNA-on-beads amplification (Positive Control: A1, A2 and A3) are higher than the Ct of cDNA in solution (C1, C2 and C3), which would be expected, given that we have already shown that the bead-based approach has a lower PCR efficiency than solution-based reactions (Fig. 3b). Evaluating the releasing efficiency of cDNA from the beads requires the corrected threshold cycle Ct s , which has been defined as previously defined. Thus, the releasing efficiency defined as (1 + E s ) −Ct /(1 + E s ) −Cts is estimated to be (mean ± s.d.) 95.6% ± 3.5% using experimentally determined cycle threshold (Ct b and Ct) and PCR efficiency (E s and E b ).

Microbeam-irradiated and bystander single cells.
Having characterized the bead-based process on chip, we were ready to integrate the PreGEM chip with the cell picker and study the single-cell response to microbeam irradiation as well as the bystander effect. 30 Individual control, 30 microbeam-irradiated and 30 bystander single cells were one by one picked from microbeam dishes via the cell picker and transferred to the microchips and processed for cell lysis, mRNA capture, RT and photocleavage. The recovered cDNAs from the on-chip reaction were then used to run digital PCR reactions with analysis on the QuantStudio platform. The QuantStudio uses limiting dilution of the reaction mix to count individual molecules, providing absolute quantification of gene expression levels using Poisson statistics. Expression levels of radiation responsive genes CDKN1A, GDF15 and PTGS2 were measured to demonstrate the utility of our single-cell preprocessing pipeline in radiation studies. CDKN1A and GDF15 are representatives of p53-regulated radiation response genes previously shown to respond predominantly in directly irradiated cells and not in bystanders 39,41 , while PTGS2 as one of the most broadly studied genes in bystander studies 42-45 is a representative of NFκB-regulated radiation www.nature.com/scientificreports www.nature.com/scientificreports/ response genes, and has been shown to respond almost identically in directly irradiated and bystander cells 39,41 . Figure 4 presents the distribution of individual control cells, bystander cells and irradiated cells based on the quantities of PTGS2, CDKN1A and GDF15 in each cell. One interesting observation is the elevated variability in expression levels in both bystander cells and irradiated cells compared to controls. For example, the expression levels of CDKN1A among irradiated cells distribute in a range of 15 to 165 counts while bystander cells have a narrower range of 0 to 45 counts, compared to that the control cells fall within a range of 0 to 25 counts. Also, it was found that the expression levels of both CDKN1A and GDF15 in irradiated cells are significantly higher than the expression levels for both genes in bystander cells (p < 0.001). The average fold-changes for CDKN1A and GDF15 in irradiated cells are 7.56 and 7.37, respectively, while the average fold-changes for these two genes in bystander cells are 2.54 and 2.19, with respect to the mean value of the two genes in control cells. While for PTGS2, the average fold-changes are 6.77 and 5.24 in irradiated cells and bystander cells, respectively. www.nature.com/scientificreports www.nature.com/scientificreports/ In order to compare the expression levels of the three genes within individual cells, regression analyses were applied to data points corresponding to expression levels of CDKN1A compared with those of PTGS2 and GDF15. We found that a linear model was not significant to explain the relationship between the expression levels of CDKN1A and PTGS2, which represent two different signaling pathways, in individual control, bystander and irradiated cells (Fig. 5a,c,e). Conversely, expression of the two p53-regulated genes, GDF15 and CDKN1A, was significantly correlated within individual cells (p < 0.05). Additionally, the r values were calculated in terms of the expression levels of CDKN1A and GDF15 for sets of control, bystander and irradiated cells. A strong linear www.nature.com/scientificreports www.nature.com/scientificreports/ relationship was found for CDKN1A against GDF15 among irradiated cells (r = 0.80), while the correlation decreased between these genes in bystander cells due to elevated variance (r = 0.65).

Discussion
To support our microbeam studies of the radiation bystander effect, we have developed a system to perform single-cell gene expression measurements following irradiation with the microbeam. This system empowers our single-cell irradiation studies, allowing individual microbeam-irradiated cells or non-hit bystander cells to be tracked, lysed, and prepared for analysis. Our approach could also be applied to studies of rare cells in many other contexts. Gene expression variability in rare cells is of great interest for applications including cancer diagnosis 46 , cancer drug development 47 and prenatal diagnosis 48 . While much can be learned through studying gene expression response in bulk cell populations, such analysis of gene expression dilutes the contribution of these rare cells to the pattern of the population. Information that covers the diversity and complexity of the cells as well as their unique molecular signatures is thereby lost.
The utility of the integrated preprocessing system was demonstrated by interfacing with QuantStudio digital PCR. The end-to-end analyses from microbeam irradiation to digital PCR gene expression measurements indicate cell-to-cell variability among individual control, microbeam-irradiated as well as bystander cells. It is critical to know that such variability is mainly due to the cell-to-cell differences in gene expression rather than being introduced by measurement noise. As can be seen from Fig. 4, the mean quantities of each gene in irradiated and bystander cells are greater than the mean of the gene in individual control cells. While the measurement noise would be expected to decrease with elevated abundance of transcripts, the variance in irradiated and bystander cells is larger than the variance in controls. This would indicate that there is minimal effect on dispersion from measurements and the variance of the single-cell data reflects actual cell-to-cell differences.
The response of lung fibroblast cells to radiation has been extensively studied by microarray and qRT-PCR 11,49-51 . For example, global gene expression four hours after bystander and direct alpha particle exposure of lung fibroblast cells was measured, revealing disparate roles of two major radiation response pathways, p53 and NFκB, in the bystander response. Activation of the p53 response pathway was found to be minimal in bystander cells, in contrast with the NFκB pathway, which appeared to respond identically in bystander and irradiated cells 39 . Additionally, lung fibroblasts were exposed to alpha-particle radiation at a dose range of 0-1.5 Gy, after which the GDF15 gene, regarded as a marker for lung injury, was assessed at the protein level and 3-fold higher expression levels were found in exposed cell culture media 24 hours after exposure 52 . The single-cell data presented in Fig. 4 are consistent with these studies, which could be explained by the previous observation that p53-regulated genes like CDKN1A and GDF15 were expressed at elevated levels in directly irradiated cultures, while showing little or no change in bystanders 39 . In contrast, the response pattern of NFκB regulated genes such as PTGS2 was found to be virtually identical in bystander and irradiated cells.
What was missed in the above studies was the heterogeneity of responses within bystander and irradiated populations. While studies in populations gave an averaged cellular response, it is clear that the presented genes were not uniformly expressed across even a small population of cells (n = 30) within each treatment set. Within irradiated cells, expression levels for any particular gene may not be uniformly elevated among the cells, in addition to inherent cell-to-cell differences revealed in the controls (Fig. 4). Signaling and response mechanisms in bystander cells are complex, involving both direct contact and gap junctions between bystanders and irradiated cells 53,54 as well as communication through extra-cellular signals in the culture media 55,56 . Thus, it is possible that the distance between an irradiated cell and a non-hit bystander, the concentration gradient of signaling molecules in the media, and the movement of the media could contribute to the observed differences in gene expression across bystander cells.
Within the same signaling pathway, CDKN1A and GDF15 both are p53-induced genes 57 , and it is clear from Fig. 5h that they are significantly correlated in microbeam-irradiated cells, and to a lesser extent in control and bystander cells. It is known that GDF15, which contains two p53 binding sites in its promoter region, is a direct target gene of p53 58 . Radiation-induced DNA damage can induce GDF15 expression in a p53-dependent manner. It has also been found that GDF15 holds a moderate but significant association with p53 after DNA damage 59 , which would indirectly reveal the correlation of genes CDKN1A and GDF15 within the p53 pathway.
Next consider the pair of CDKN1A and PTGS2, which represents activation of different pathways including p53 and NFκB. The former gene plays essential roles in the DNA damage response 60 and the latter counteracts p53 activity as well as inhibits DNA damage-induced apoptosis 61,62 . As seen from Fig. 5, the correlations between the two genes in controls and bystanders are weak and positive while the two genes were found to be negatively related in irradiated cells. The negative trend observed between CDKN1A and PTGS2 in irradiated cells may imply a transcriptional cross talk between NF-κB and p53, in which both pathways inhibit each other's ability to stimulate gene expression and this process is likely dependent on the relative levels of the two transcripts 63 . It has been suggested that stimulation of NF-κB promotes resistance to programmed cell death while the activation of p53 is associated with the induction of apoptosis or cell cycle arrest following DNA damage 64,65 . If this is indeed the case, it is possible that NF-κB activation occurs simultaneously with the induction of p53 in a few individual cells shown in Fig. 5e, which exhibited equally high levels of expression of both genes, but that activation of the two pathways is not tightly linked within the same cell.

conclusions
The data presented here demonstrate the integrated preprocessing approach for single-cell irradiation studies, coupled to qPCR or digital PCR, provides a powerful method to investigate the variability of gene expressions in rare cells with microbeam irradiation, and many other applications where rare cells are of interest.  20 and allowed to attach for at least 1 h (Fig. 6). Immediately before irradiation, the culture medium was removed from the microbeam dish. The microbeam dish was positioned adjacent to the microbeam exit window and a moisture containment collar was placed over the objective lens to keep the cells from becoming dehydrated during the approximately 15 min irradiation time. The RARAF microbeam system has integrated image analysis, which was used to visualize the nuclei stained with Hoechst 33342. A computer automated nuclear irradiation protocol found the location of every stained nucleus on the dish. We used the RARAF 5 MV Singleton Accelerator to deliver five doubly charged helium particles (He 2+ ), simulating an alpha particle, with an energy of 6.0 MeV to the center of each stained nucleus. The ion beam was measured to have a width of <5 µm for all irradiations, and based on SRIM calculations (Stopping Range of Ions in Matter) 66 , each He 2+ particle had a linear energy transfer (LET) of approximately 100 keV/µm. Cells showing nuclear stain were microbeam-irradiated while the MitoTracker Green-stained cells were not targeted and therefore were designated as "bystander" (Fig. 6). The particle fluence was measured by a gas proportional counter positioned above the cells. After every nuclear stained cell on the plate had been irradiated the dish was removed from the irradiation endstation, fresh medium was added to the dish, and the dish was placed in the incubator at 37 °C until ready for cell picking at 4 hours post irradiation. All control cells were stained with Hoechst 33342 and sham irradiated.
Cell picking. Cells were viewed on a custom microscope using a 10x objective with live imaging (pco.edge 5.5, PCO, Kelheim, Lower Bavaria, Germany). Custom software allowed for imaging and lighting control as well as computer adjustment of the position of the pipette mounted on a 6 axis stage (X, Y, Z, azimuth, elevation and retraction; Zaber Technologies Inc., Vancouver, BC). 40 uL of 0.25% trypsin-EDTA was added to the microbeam dish prior to picking to aid in cell release. To pick a cell, the pipette was positioned directly adjacent to a cell of interest on the microbeam dish. Multi-color illumination, with the use of a Lumencor Aura Light Engine (Lumencor, Beaverton, OR), a light source (Morrell Instruments, Mellville, NY) and a camera (PCO-Tech Inc, Romulus, MI), allowed for simultaneous imaging of both the nuclear stain and cytoplasmic stain, thus permitting determination of the cell as "irradiated" or "bystander", respectively. Single cells were aspirated into the micropipette along with 1-3 µL of 0.25% trypsin-EDTA. After picking a cell, the micropipette was moved out of the microbeam dish and positioned within the input well of a PreGEM chip, where the picked cell was ejected directly into lysis buffer (Thermo Fisher, Waltham, MA). www.nature.com/scientificreports www.nature.com/scientificreports/ Fabrication of PreGEM chip. The chip (Fig. 7a) is fabricated using standard soft lithography microfabrication techniques 67 . SU-8 photoresist (MicroChem Corp., Newton, MA) was spun-coated on a silicon wafer and patterned with photolithography to define the liquid channel mold. Then, PDMS (Dow Corning) was poured over the mold and an additional vapor barrier 33 was embedded in the PDMS. Sheets bearing the microfluidic features were then peeled off the mold followed by inlet and outlet hole punching. The molded PDMS was sealed to a glass slide using oxygen plasma bonding.
On-chip RNA capture. Different copy numbers of XenoRNA template (2 × 10 4 , 5 × 10 4 and 10 5 , Thermo Fisher Scientific Inc.) were introduced to the microchip, and were captured by 3.5 × 10 6 beads (Thermo Fisher Scientific Inc). Separating the magnetic beads from the remaining solution using a magnet, the effluent solution was transferred to micro tubes and mixed with another solution containing 3.5 × 10 6 beads, which was followed by off-chip RT. In the end, 40-cycle qPCRs were run to detect the amount of un-captured templates.
On-chip RT characterization. The efficiencies of the on-chip bead-based RT were investigated using various copy numbers of XenoRNA templates (2 × 10 4 , 5 × 10 4 and 10 5 , Thermo Fisher Scientific Inc.). The on-chip RTs were first conducted after mixing RT reagents, XenoRNA, biotinylated oligo(dT) 25 primer (Integrated DNA Technologies Inc.) and 3.5 × 10 6 streptavidin magnetic beads. Subsequently, PCRs were performed with the products of RT using the qPCR system. In comparison for the on-chip efficiency, in-tube RTs were performed using a thermocycler followed by qPCR characterization. Additionally, to characterize the effects of beads on reaction efficiency, bead-based PCRs were subsequently performed with the products of in-tube RT by mixing and binding biotinylated cDNA to streptavidin beads, comparing the bead-based approach with the standard bench-top solution-based approach.
Single-cell processing on the PreGEM chip. The input well of the PreGEM chip was pre-filled with the lysis buffer, so as to immediately lyse the input cell after it was picked and placed within the input well. Next, mRNA templates in the cell lysate were collected by photocleavable 5′-PC Biotin-(dT) 25 -3′ bead tether, which was made by incubating magnetic streptavidin beads (Thermo Fisher Scientific Inc.) with photocleavable biotinylated oligo(dT) 25 (Integrated DNA Technologies Inc.). The principle of mRNA capture relies on base pairing between the polyA tails of the mRNA and oligo(dT) 25 immobilized on the surface of the beads. In the reaction of reverse transcription (RT), the bead-bound oligo(dT) 25 functions as a primer for synthesis of cDNA. Extracted mRNA were moved from inlet well to reaction chamber using a magnet, followed by the introduction of reverse transcription mix (dNTP, RT buffer, MgCl 2 solution, Reverse Transcriptase, RNase inhibitor). After RT (10 min at 25 °C and 50 min at 42 °C), the microchip was exposed under a UVA lamp (Sylvania Inc) photocleaving the PC linker, which released the cDNA from the beads. In-tube RT, qPCR and digital PCR (dPCR). In-tube RTs were run using a thermocycler (Eppendorf) with a standard protocol (10 min at 25 °C and 50 min at 42 °C). qPCRs were performed using 96-well reaction plates run on the 7900HT Fast Real-Time PCR (Applied Biosystems). 5 μL RT products, 6.25 μL nuclease-free water, 1.25 μL Taqman gene expression assay (Hs99999905_m1, Thermo Fisher Scientific Inc.) and 12.5 μL Taqman gene expression master mix were loaded on the plates, which were run for 45 cycles. The QuantStudio 3D digital PCR system (Thermo Fisher Scientific Inc.) was used to conduct digital PCRs. Starting samples in a final volume of 15 μL that consisted of bead-excluded RT product, Taqman gene expression master mix (Thermo Fisher Scientific Inc.), Taqman gene expression assay (Hs99999142_m1, Hs00153133_m1, Hs00171132_m1, Thermo Fisher Scientific Inc.) and nuclease-free water were loaded onto Digital PCR Chips (Applied Biosystems) via the QuantStudio 3D Digital Chip Loader (Applied Biosystems). The chips were immersed in oil to avoid evaporation and sealed with UV sealant (Applied Biosystems). This was followed by running 40 cycles of dPCR using GeneAmp ® PCR system 9700 Thermocycler, which is adapted for digital PCR reactions. The chips were then loaded onto the QuantStudio ™ 3D Digital PCR Instrument, and the end-point fluorescence of each well on the chips was measured using QuantStudio ™ 3D Digital PCR software v3.0. The fluorescence data were read and analyzed using QuantStudio 3D Analysis Suite Cloud Software.
Statistical analysis. Statistical differences between two groups of cells were determined using unpaired two-sample t-test or one-way ANOVA followed by post-hoc Tukey test, where necessary. Pearson correlation coefficient (r value) was used to measure the linear relationship between expression levels of two genes. Defined as the covariance of two variables divided by the product of their standard deviations and calculated as follows: , where x and y represent variables and n is the number of experiments. The Kolmogorov-Smirnov test was used to test the normality of the samples and to decide whether the distribution functions of two samples are significantly different. R (version 3.1.1) was used to perform statistical analysis.