1. Introduction

Developmentally appropriate genetic signaling depends on internal feedback networks responding to external signals. Recent works have demonstrated that single-cell measurements are key to understanding how these networks respond to external factors because single-cell behavior is rarely normally distributed about a mean value ^[1,2]. This is partially due to RNA bursting, during which RNA are stochastically produced in the nucleus and after a series of steps arrive in the cytoplasm for translation ^[1,2]. This spatial separation provides information on the timing of events, even for populations of fixed cells. For example, Neuert et al. measured the single-cell RNA expression of thousands of budding yeasts in response to osmotic shock for multiple genes. Using this dense dataset, they created in silico models that predicted how the internal regulatory networks would respond to new external conditions or genetic perturbations to the network itself ^[3]. This insight would not be possible using standard ensemble methods because the variability in spatial distribution of RNA and copy number variation in RNA would be lost due to population averaging ^[1,2,3,4]. This set of techniques provides a powerful tool to image fixed cell populations and infer the temporal sequence of events within a signaling network.

Vascular development is partially driven by highly proliferative endothelial progenitor cells (EPC). Numerous assays have been utilized to enumerate and isolate EPC by flow cytometry and primary cell culture ^[5,6,7]. Although conflicting definitions exist in the literature, Critser and Yoder concluded in 2010 that late outgrowth endothelial colony-forming cells (ECFC), which demonstrate the ability for self-renewal and neo-angiogenesis, both in vitro and in animal models, can uniquely be described as progenitors of endothelial cells ^[6,8]. A key feature of ECFC is tube formation when plated on Matrigel or cross-linked collagen, allowing for the study of in vitro angiogenesis without the need for complicated growth matrices. We previously demonstrated that ECFC are reduced in prematurely born infants that go on to develop bronchopulmonary dysplasia (BPD) and infants born to mothers with preeclampsia ^[9,10,11]. For prematurely born infants the introduction to room air after birth leads to an increase in environmental oxygen at a developmentally abnormal time.

One potential regulator of ECFC environmental response to hyperoxia is Thymosin β4 (Tβ4). Cytosolic Tβ4 protein mediates the response of endothelial cells to changes in oxygen environment by modulating the stability of Hypoxia Inducible Factor 1α (HIF-1α). HIF-1α is historically associated with cellular response to hypoxia, not hyperoxia. Recent work has demonstrated a role for HIF-1α in signaling preservation across all changes in oxygen environment, particularly in endothelial cells ^{[12,13,14,15,16]}. One key feature of this signaling network is a self-modulating feedback mechanism for Tβ4. Jo et al. demonstrated altered HIF-1α nuclear localization as a function of increasing Tβ4 expression in endothelial cells ^[12]. Based on these results, we applied high-throughput imaging to determine if changes in the spatial localization of Tβ4 and HIF-1α in ECFC are associated with changes in RNA copy number of vascular endothelial growth factor (vegf) and endothelial nitric oxide synthase (eNOS) at the single-cell level.

vegf and eNOS are potent pro-proliferation and pro-angiogenesis genes that are associated with many developmental diseases in prematurely born infants ^{[9,17,18,19,20]}. Our laboratory previously showed that ECFC derived from the umbilical cord blood of prematurely born infants (pt-ECFC) demonstrated markedly different proliferation in vitro as compared to ECFC derived from the umbilical cord blood of full term birth infants (t-ECFC) ^[9,10,11]. Based on these results and others, we hypothesized that exposure to hyperoxia introduces a specific environmental signal that blunts pro-proliferation and pro-angiogenesis genetic signaling networks in a sub-population of pt-ECFC. Specifically, we hypothesized that elevated basal cytosolic Tβ4 protein concentration in a subpopulation of pt-ECFC is associated with diminished HIF-1α nuclear localization and diminished early-time (within the first 60 minutes) vegf and eNOS RNA production in response to hyperoxia as compared to t-ECFC.

To test this hypothesis, we utilized high-throughput single-cell imaging of protein and RNA expression to quantify pt-ECFC and t-ECFC grown in vitro room air and hyperoxia. We quantified Tβ4 and HIF-1α protein expression using immunofluorescence as well as vegf and eNOS RNA expression using single-molecule fluorescence in-situ hybridization (smFISH ^[21,22]). Because these techniques require cell fixation we were unable to track single-cell temporal correlations. We instead built statistical spatial measures of how independent populations of pt-ECC or t-ECFC respond to either room air or hyperoxia treatment by measuring cell populations across multiple time points and multiple environmental conditions. This allowed us to infer potential casual relationships through statistical analysis.

We found that pt-ECFC have elevated basal cytosolic Tβ4 as compared to t-ECFC. This was statistically associated with lower HIF-1α nuclear localization. During the first hour of in vitro hyperoxia growth, distinct single-cell populations emerged in pt-ECFC that are not present in t-ECFC. Based on these findings we speculate that this pathway may directly control the pro-angiogenic response of ECFC and propose a set of future experiments.

---

2. Materials and Methods

---

2.1. Cord blood collection and ECFC isolation

Cord blood (CB) collection and endothelial colony forming cells (ECFC) were prepared as previously described ^[23]. CB samples were obtained following informed consent from term (gestational age greater than 37 weeks) and preterm (gestational age less than 36 weeks) infants at the University of Colorado Anschutz Inpatient Pavilion. The Colorado Multiple Institutional Review Board approved all protocols. All CB was processed within 24 hours of collection.

Mononuclear cells (MNCs) were plated on cell culture plates coated with type 1 rat-tail collagen (BD Biosciences; San Jose, CA) at a density of 5 × 10⁶ cells/cm². Complete EGM-2 medium (Lonza, CC-3612; with 10% fetal bovine serum) was changed daily for seven days. After seven days, media were changed three times per week and ECFC colonies were enumerated on day 14 using light microscopy. To ensure that ECFC were not migrating or being removed during feeding, the removed media from 6 cell wells was re-plated on collagen-coated plates and observed for 14 days. Cells were then collected into cryogenic tubes and frozen under liquid nitrogen until future studies were performed.

---

2.2. ECFC in vitro culture

ECFC were removed from cryo-suspension and studied at passage 3–5 in complete EGM-2 media. Individual #1 cover slips (Corning, 2845-25) were cleaned in ethanol followed by 1 × PBS. ECFC were plated onto coverslips at low density and cultured in either room air with 5% CO₂ or 50% O₂ with 5% CO₂. At 0, 15, 30 or 60 minutes the EGM-2 medium was replaced with 3.7% paraformaldehyde (Sigma, 47608) for 15 minutes, followed by 1 × PBS (Sigma, P5493). All steps were performed within the incubator to ensure a near constant oxygen environment. Following completion of all coverslips, the PBS in each well was exchanged for 70% ethanol (Sigma, 277649).

---

2.3. High-throughput single-molecule, single-cell imaging

Automated imaging was performed using a custom-built microscope, consisting of an Olympus IX71 microscope base, X-Y translation stage (Mad City Labs Microstage), objective piezo (Mad City Labs F200S), oil-immersion 100 × NA1.3 objective (Olympus UPLSAPO 100XO), multi-color LED light source (Lumencor Spectra-X), LED light source specific filter set (Semrock LED-DA/FI/TR/Cy5-4X-A), and a sCMOS camera (Hamamatsu C11400-22CU). Automation and acquisition were performed using MicroManager v1.4 running on a Windows 7 64-bit laptop ^[24]. Analysis was performed on a computing cluster with 48 cores, 128 gigabytes of RAM, 2 NVidia TITAN GPUs, and hot swappable storage (Titanus Computers). Custom analysis software was used to analyze all data in both ImageJ ^[19] and MATLAB (Mathworks).

Fifty image areas, containing at minimum two cells per area, were identified using the automated stage and objective piezo. At each image area, a simple contrast based autofocus routine implemented in MicroManager was used to correct mechanical or thermal drift. Forty axial image planes were acquired at 250 nm spacing, centered around the algorithm selected best focus point. At each focal plane, four independent images were captured, one each for DAPI, Alexa 488, Alexa 561, and Alexa 647 fluorescent labels.

---

2.4. Tβ4 and HIF-1α fluorescent labeling

All fluorescent labeling steps were adapted from previous work by our group and others ^[20,21,22]. Following overnight permeabilization in 70% ethanol at 4 ℃, coverslips were washed in 2 mL of 1 × PBS at room temperature followed by incubation in blocking buffer (Thermo-Fisher, 37515). Coverslips were incubated in primary antibody (1:100 dilution in 1 × PBS) for Tβ4 (Abcam, ab14335) and HIF-1α (Abcam, ab51608) overnight, washed in 2 mL of 1 × PBS, incubated in secondary antibody (1:100 dilution in 1 × PBS; Abcam, ab175471; Abcam, ab150115) at 37 ℃, and washed in 2 mL of 1 × PBS with 1 drop of pre-mixed DAPI (Life Technologies, R37606).

The following steps were performed in low light to prevent bleaching. Coverslips were hybridized with 20 μL of GLOX and Phalloidin-Alexa488 (Life Technologies, A22287) for 30 minutes at room temperature and then washed with GLOX buffer for 15 minutes at room temperature. Coverslips were plated on standard microscope slides with 8 μL GLOX plus enzyme (imaging buffer) and nail polish was applied around the coverslip to seal.

---

2.5. Tβ4 and HIF-1α image analysis

---

2.5.1. GPU-based de-convolution

Image de-convolution was performed in ImageJ, utilizing software provided by the Butte Group ^[25], prior to any further image processing steps.

---

2.5.2. Cell and nuclei identification

Custom in-house MATLAB software was utilized to filter and identify individual cells in every image prior to protein analysis ^[4,26,27]. The user was presented a maximum intensity projection image of the nuclei and actin cytoskeleton fluorescent labels. Each cell was manually circled and recorded until all cells in an image are identified. A watershed algorithm was used to automatically detect the nucleus within each identified cell. Overlapping cells are not included in further data analysis to avoid over-counting. Cell and nuclei outlines are saved as text-files for use later in the analysis.

---

2.5.3. Protein detection and quantification

The basic workflow of protein detection and quantification was: 1. image filtering and 2. corrected total cell fluorescence (CTCF) calculation. In-house MATLAB and ImageJ software were utilized to median filter and then calculate the CTCF value for every cell. The background for CTCF was independently determined for each image.

---

2.6. vegf and eNOS RNA labeling

---

2.6.1. Single-molecule fluorescence in-situ hybridization

All hybridization steps were adapted from previous work by our group and others ^[4,21,22]. Following the overnight permeabilization, coverslips were washed in 2mL of 10% wash buffer for 30 minutes at 37 ℃. Coverslips were incubated with 20 μL of vegf and eNOS probe solution (1:1:1000 in 10% hybridization buffer) for 8 hours at 37 ℃. Both vegf and eNOS probe sets were designed using both in-house and commercially available probe designing technologies (Biosearch Technologies). Following hybridization, coverslips were washed with 10% wash buffer for 30 minutes at 37 ℃ and then incubated in 10% wash buffer with 1 drop pre-mixed DAPI solution (Life Technologies, R37606) in each well for 30 minutes at 37 ℃. Coverslips were then washed with GLOX buffer for 15 minutes at room temperature.

---

2.6.2. Hybridization buffers

10% wash buffer: 5 mL deionized formamide (Ambion, AM9342), 5 mL 20 × SSC (Ambion, AM9763), 40 mL nuclease-free water (Ambion, AM9932). Hybridization buffer: 1 mL deionized formamide (Ambion, AM9342), 10 mg E. coli tRNA (Roche, 10109541001), 1 mL 20× SSC (Ambion, AM9763), 40 μL ultrapure BSA (Ambion, AM2618), 1 g dextran sulfate (Sigma, D8906-50G), 100 μL of 200mM vanadyl-ribonucleoside complex (Sigma, 94742-1ML), nuclease-free water (Ambion, AM9932) to 10mL total solution. GLOX buffer: 100 μL 1M Tris (Ambion, AM9855G), 1mL 20× SSC (Ambion, AM9763), 400 μL of 10% Glucose, 8.5mL nuclease-free water (Ambion, AM9932).

---

2.6.3. Imaging preparation

The following steps were performed in low light to prevent bleaching. Coverslips were hybridized with 20 μL of GLOX and Phalloidin-Alexa488 (Life Technologies, A22287) for 30 minutes at room temperature and then washed with GLOX buffer for 15 minutes at room temperature. Coverslips were plated on standard microscope slides with 8 μL GLOX plus enzyme (imaging buffer) and nail polish was applied around the coverslip to seal.

---

2.6.4. Imaging buffer

GLOX plus enzyme: 1 μL 3.7 mg/mL Glucose Oxidase (Sigma, G7141), 1 μL Catalase (Sigma, C30), 100 μL GLOX buffer.

---

2.7. vegf and eNOS image analysis

---

2.7.1. GPU-based de-convolution

Image de-convolution was performed in ImageJ, utilizing software provided by the Butte Group ^[25], prior to any further image processing steps.

---

2.7.2. Cell and nuclei identification

FISH-QUANT and custom in-house MATLAB software were utilized to filter and identify individual cells in every image prior to RNA analysis ^[4,26,27] The user was presented a maximum intensity projection image of the nuclei and actin cytoskeleton fluorescent labels. Each cell was manually circled and recorded until all cells in an image are identified. A watershed algorithm was used to automatically detect the nucleus within each identified cell. Overlapping cells were not included in further data analysis to avoid over-counting. Cell and nuclei outlines are saved as text-files for use later in the analysis.

---

2.7.3. RNA detection and quantification

The basic workflow of RNA detection and quantification was: (1) image filtering, (2) spot pre-detection, (3) spot fitting, (4) non-specific versus specific determination, (5) spot assignment to individual cells, (6) batch processing.

FISH-QUANT was utilized to identify, fit, and quantify individual RNA labeled using smFISH. The details of FISH-QUANT are described in a prior publication ^[26]. Background fluorescence is typically due to FISH probes that were either not washed out or non-specifically bound during the labeling procedure ^[4]. FISH-QUANT utilizes a three-dimensional Gaussian to fit each possible spot and the user must discriminate background from specific binding by manual thresholding of several extracted fit parameters (e.g. lateral fit uncertainty, axial fit uncertainty, amplitude of Gaussian).

---

2.8. Statistical analysis of imaging results

Single-cell protein CTCF and RNA distributions were loaded in MATLAB for statistical analysis. Before further analysis, all measurements were divided by the maximum intensity projection of the measurement's corresponding cell compartment area (nucleus, cytoplasm, or total cell). Area normalized marginal distributions were tested for independence using the two sample Kolmogorov Smirnov (K-S) test implemented in the Statistics Toolbox for MATLAB (null hypothesis rejection at α = 0.05). Area normalized marginal distribution independence were tested at each measured time for a given protein or RNA species for all cell types and growth conditions (population response). Additionally, we tested within each protein or RNA species for a specific cell type and growth condition against measure time points (temporal response). Area normalized multivariate joint distributions were tested for independence using the non-parametric and distribution free multivariate distribution distance tests, implemented in the highdim package for MATLAB (null hypothesis rejection at α = 0.05) ^[28]. This allowed for direct comparison of the non-Gaussian distributions common in single-cell gene expression imaging experiments using a non-parametric test.

---

3. Results

---

3.1. Cytosolic basal Tβ4 expression in ECFC

To begin, we tested if pt-ECFC Tβ4 expression was elevated in multiple primary ECFC cultures. To determine Tβ4 expression, we thawed and prepared individual populations of ECFC from four term and four preterm patients. After three to five passages, we fixed each population and labeled for Tβ4 using immunofluorescence. We assessed protein content using corrected total cell fluorescence (CTCF) for the nucleus, cytoplasm, and total cell. We normalized CTCF values by the maximum projection area of the feature of interest to limit copy number variability due to cell size.

Using Tβ4 expression from at least 200 cells per patient we find that cytosolic Tβ4 protein expression is elevated across pt-ECFC (n = 4 patients) as compared to t-ECFC (n = 4 patients) at basal conditions (Figure 1). We found that total cell and nuclear Tβ4 expression are not statistically different.

---

3.2. Tβ4 protein expression and HIF-1α protein nuclear localization

After establishing that pt-ECFC have elevated Tβ4, we moved on to test if pt-ECFC display diminished HIF-1α nuclear localization. While nuclear localization is not a perfect indicator of transcription factor activity, it is a requirement for a transcription factor to be in the nucleus to activate a gene network. To determine protein localization, we cultured a mixture population of the four patients used in Figure 1 in room air or hyperoxia (50% O₂). At 0, 15, and 30 minutes of growth we fixed and labeled for both Tβ4 and HIF-1α using immunofluorescence. We assessed protein content in the nucleus, cytoplasm, and total cell using CTCF. We normalized CTCF values by the maximum projection area of the feature of interest to limit copy number variability due to cell size.

We found that HIF-1α nuclear localization was diminished for pt-ECFC grown in vitro room air or hyperoxia (top 2 rows; blue line; Figure 2A). In contrast, we found that HIF-1α nuclear localization was immediate and sustained for t-ECFC grown in vitro room air and hyperoxia (bottom 2 rows; blue line; Figure 2A). We found that there was minimal difference in Tβ4 protein localization for pt-ECFC grown in vitro room air or hyperoxia versus t-ECFC grown in vitro room air or hyperoxia respectively (top 2 rows versus bottom 2 rows; Figure 2B). This suggests that it was the initial increased concentration in Tβ4, or other factors not measured here, that lead to the observed reduction in nuclear HIF-1α localization for pt-ECFC grown in vitro hyperoxia ^[12].

---

3.3. vegf and eNOS RNA expression

One consequence of reduced HIF-1α nuclear localization was a potential reduction in vegf and eNOS RNA production. We utilized the same experimental design as above but instead labeled for RNA transcripts at each time point. We found the pt-ECFC do not begin production of vegf or eNOS RNA in either room air or hyperoxia until 60 minutes (top 2 rows; blue line; Figure 3A and B). In contrast, we found t-ECFC began and sustained vegf and eNOS RNA production almost immediately (bottom 2 rows; blue line; Figure 3A and 3B). This production resulted in increased mature RNA in the cytoplasm by 60 minutes (bottom 2 rows; red line; Figure 3A and 3B).

---

3.4. Protein or RNA expression as a function of nucleus to cytoplasm area ratio

Based on our above findings that there were differences in Tβ4, HIF-1α, vegf, and eNOS marginal distributions both between pt-ECFC and t-ECFC as well as in vitro room air and hyperoxia growth, we utilized our single-cell dataset to construct multi-dimensional distributions. Based on previous findings that the ratio of nuclear area to cytoplasmic area (N:C ratio) may be a marker of proliferating ECFC ^[29], we plotted the N:C ratio versus Tβ4 CTCF and HIF-1α CTCF (Figure 4A) and N:C ratio versus vegf RNA transcripts and eNOS RNA transcripts (Figure 4B). Unlike the marginal distributions presented in Figures 2–3, we utilized the total amount of protein or RNA in each cell to simplify the analysis.

This analysis provides a visualization of how the four different cell populations evolved over time, fully leveraging the information within this single-cell dataset. pt-ECFC grown in vitro HO (green) demonstrated different behavior from all other populations for both protein (Figure 4A) and RNA (Figure 4B) (p < 0.01, multi-dimensional distance test). A large heterogeneity in both cell size and protein or RNA expression was also apparent using this analysis, supporting the idea that ensemble measurements may be averaging over several different single cell behaviors, masking these differences.

---

4. Discussion

Our analyses suggested that previous observations of blunted proliferation of pt-ECFC in hyperoxia ^[9,23] may be due to developmentally appropriate levels of pro-growth proteins, such as Tβ4. These levels were elevated as compared to t-ECFC, most likely due to internal and external factors driving ECFC proliferation. Paradoxically, the release of Tβ4 into the cytosol in response to hyperoxia reduced HIF-1α stability in this high-Tβ4 sub-population ^[12]. That in turn lead to lowered vegf and eNOS RNA production and the emergence of a less progenitor-like population. The overall results of this network motif matched previous ensemble proliferation measurements, but the details are only observable using single-cell methodologies.

Because these data were "stop-motion" movies of distinct cell populations, we cannot confirm that those cells with elevated Tβ4 were the same cells that emerged with reduced HIF-1α nuclear localization and subsequent lower vegf and eNOS RNA expression. Future experiments that rely on live-cell fluorescent reporters of Tβ4 and HIF-1α protein followed by smFISH measurement of vegf and eNOS RNA may provide key insight into the temporal evolution of these heterogeneous cell populations.

---

5. Conclusion

In this study, we attempted to address our hypothesis that elevated basal cytosolic Tβ4 protein concentration in a subpopulation of pt-ECFC diminishes early-time (within the first 60 minutes) vegf and eNOS RNA production in response to hyperoxia as compared to t-ECFC. Using single-cell imaging, we found statistical evidence that elevated cytosolic Tβ4 in pt-ECFC were associated with diminished HIF-1α nuclear localization and diminished vegf and eNOS RNA copy number. Our data and analyses, combined with previous studies, suggest that Tβ4 may play a key role in determining ECFC vegf and eNOS response to changes in oxygen environment. Our data additionally suggest a role for HIF-1α, further expanding the scope of this key transcription factor to the ECFC hyperoxia response.

---

Acknowledgements

M.A.A. and D.P.S acknowledge startup funding from the University of Colorado Denver College of Liberal Arts and Sciences. G.J.S. acknowledges funding from the National Institute of Health (NHLBI HL68702). C.D.B. acknowledges funding from the National Institutes of Health (K23 HL121090-01A1).

---

Conflict of Interest

All authors declare no conflicts of interest.

---