Image3C is an imaging tool developed to study tissue composition at single-cell resolution in research organisms for which antibodies and pre-existing knowledge about cell types are not readily available (Figure 1, Table 1). Image3C allows for high-throughput and unbiased analysis in scenarios where manual counting and observer-based cell identification are currently the only options. Image3C includes all the components required for compensating captured images, quantifying multiple features for each event, clustering the events, visualizing and exploring the data, and training and using the CNN for analyzing subsequent samples and integrating multiple experiments (Figure 1, Figure 1—figure supplement 1).

Once a single cell suspension is prepared from the organism of interest, the cells are stained with a combination of dyes that are expected to function irrespectively of the species used and which have high affinity for specific cellular organelles such as nuclei or molecules associated with metabolic states such as reactive oxygen species (ROS). We validated reagents experimentally by determining that nuclear dyes stain intracellular material matching expected characteristics of nuclear DNA or by activation of cells with drugs to change their metabolic state. The labeled samples are then run on the ImageStream^X Mark II (Figure 1A). ImageStream is a commercially available image-based flow cytometer, whose diffusion in laboratory settings is increasing and that provides highly reproducible images of cells that can be compared across days of acquisitions and experiments. For this approach, no microfluidic devices or custom-made and highly specialized microscopes are required (Table 1), and, if desired, the users can test the Image3C pipeline also on images acquired at standard microscopes, remembering to carefully control for batch effects.

Once images of individual events are collected for each channel of interest, feature values from both morphological and fluorescent data, such as cell size and nuclear size, are extracted from the cell images using IDEAS software (Amnis Millipore, free for download upon creation of Amnis user account) (Figure 1B; see Supplementary files 1 and 2 for feature description). Correlation between features is calculated and redundant features are trimmed as well as samples that, among replicates, are outliers (Figure 2—figure supplements 1 and 2). This prevents clustering artifacts potentially caused by having multiple features providing the same information or including samples that are not representative (Figure 1B). During this step, while the number of features was usually reduced significantly, the correlation between replicates was always high and outliers were rarely observed (Figure 2—figure supplements 1 and 2). Finally, fluorescence intensity features are transformed to improve homogeneity of variance of distributions and, if used, DNA staining is normalized to remove intensity drift between samples and thus align the 2N and 4N DNA content histogram peaks (Figure 1B, Figure 2—figure supplement 3).

Exported feature quantifications are used for clustering the events. Dimensionality reduction and visualization of clusters is achieved by generating force-directed layout (FDL) graphs in the VorteX clustering environment (Figure 1C) (free to install) (Samusik et al., 2016). Cell images for events within each cluster can be visualized using FCS Express Plus together with custom R scripts (Figure 1D). These visualization tools and the cluster feature averages (i.e., the mean value of each feature for each cluster) (Figure 1E) allow to explore the images of selected groups of events and the features that differ between cells belonging to separate clusters. If control and treatment samples are included, a statistical analysis using negative binomial regression to compare cell counts per cluster between samples is also available in the Image3C pipeline. This high-throughput and unbiased analysis provides a comprehensive understanding of a cell population composition at higher resolution than what is possible with traditional histological methods.

Once this pipeline is run on a first set of samples (e.g., homeostatic state) and the cell clusters are defined for the tissue of interest, the images and the relative clustering IDs can be used to train a CNN classifier in an unbiased way (Figure 1F), including the ability to score frequency of ‘new’ cell types that do not match any of the clusters identified at homeostasis. Therefore, future experiments in the same tissue used for training the CNN classifier can be analyzed directly through the CNN (Figure 1G). This significantly reduces the number of steps and time required to process data collected subsequently. An even greater advantage is represented by the fact that, in the absence of CNN, every time new experimental sets are run it would be necessary to go again through the de novo clustering part of the pipeline (Figure 1B–E) and the new set of clusters would need to be cross-annotated to be compared with cell population composition observed in previous experiments. Manually matching clusters between different experimental sets might be a source of errors, mainly if the user is not familiar with the cell types present in the sample and if specific biomarkers or pre-existing knowledge about cell types and morphology are not available. The CNN splits all the cell images in the classes defined during the training step and allows to compare the abundance of cells with same morphology among different samples without the need to cross-annotate clusters (Figure 1F, G). The CNN inclusion in Image3C and the reproducibility of image acquisition through image-based flow cytometry allows use of the clusters defined from one experiment (e.g., homeostatic state) to set up a classifier in an unbiased way for later use as a reference in analyzing the effects of experimental manipulations on these cell populations. We conclude from these results that Image3C can perform de novo high-throughput characterizations of population composition and define specific cell type changes between homeostatic and experimentally perturbed samples across multiple experiments.

To test whether Image3C could identify homogeneous and biologically meaningful cell populations, we used the research organism Danio rerio. We obtained cells from adult female zebrafish WKM (location of the hematopoietic tissue) in homeostasis condition, stained them, and ran on the ImageStream^X Mark II. We analyzed intrinsic morphological and fluorescent features, such as cell and nuclear size, shape, and darkfield signal (side scatter [SSC]). Feature values were extracted from each cell image and processed through our pipeline (see Supplementary file 1 for feature description). Clustering by the final set of normalized and non-redundant morphological and fluorescent features produced distinct cell populations (Figure 2A–C, Figure 2—figure supplements 1–3).

Figure 2 with 3 supplements see all

Download asset Open asset

Image3C can distinguish between the major classes of cells present in zebrafish WKM (Figure 2; Supplementary files 3 and 4) that were described using standard sorting flow cytometry and morphological staining approaches (Traver et al., 2003). It is noteworthy that Image3C can clearly identify dead cells and debris (Figure 2A, B) allowing to optimize experimental protocols in order to minimize cell death and run the subsequent analysis only on the intact, live cells. Image3C can identify cells with outstanding morphological features, such as neutrophils from other myelomonocytes (Figure 2B, C). Based on zebrafish neutrophil characteristics such as high granularity, high complexity, and low circularity of the nuclei (Lugo-Villarino et al., 2010), this type of granulocytes can be easily distinguished. Other types of myelomonocytes, such as monocytes and eosinophils, are here merged in the same cluster since in zebrafish they share many morphological characteristics (Lugo-Villarino et al., 2010). Similarly, using only intrinsic morphological features for the clustering, different lymphocytes (B and T-cells) and hematopoietic stem cells cannot be separated from each other, but they can be clearly distinguished from the myelomonocytes (Figure 2A, B). Within the lymphocytes/progenitors fraction, we find two clusters (Dr1 and Dr7) that mainly differ in cell diameter (Figure 2C). Whether this morphological difference has a biological implication needs to be studied in future experiments.

Image3C also enables the quantification of cell populations (clusters or CNN classes)-relative abundance, an important tool for comparing population composition across different treatment groups under different environmental conditions (Peuß et al., 2020). Here, we compared our results with previously published data to validate our method. Although a direct comparison with results from classical approaches (Traver et al., 2003) is not possible since we gated out (removed analytically) mature erythrocytes before clustering (Materials and methods), the myelomonocyte to lymphocyte ratio (M/L ratio = 1.59) is similar to the one obtained with classic histological approaches (mean M/L ratio = 1.35) (Figure 2D; Traver et al., 2003).

Next, we sought to determine whether Image3C could be used to characterize and quantify biological processes by identifying a tissue of interest and then comparing cellular composition dynamics, functions, and physiological responses of specific cell types across a range of experimental conditions. Our goal was to detect statistically significant changes in cluster relative abundance between control and treated samples to gain a more detailed understanding of cell population dynamics and individual cell function.

As proof of concept, we performed a standard phagocytosis assay using WKM tissue from female adult zebrafish. The single cell suspension was incubated with CellTrace Violet labeled Staphylococcus aureus (CTV-S. aureus) and with dihydrorhodamine-123 (DHR), a ROS indicator that becomes fluorescent if oxidized to report oxidative bursting following phagocytosis. As controls, we inhibited phagocytosis through cytoskeletal impairment with cytochalasin B (CCB) incubation or through incubation with bacteria at lowered temperature by placing culture plates on ice. Events collected on the ImageStream^X Mark II were analyzed with Image3C and clustered in 26 distinct clusters using quantifications of morphological and fluorescent features (see Supplementary file 2 for feature description), including nuclear staining, phagocytized S. aureus, and DHR positivity (Figure 3A, Figure 3—figure supplement 1). Professional phagocytes are defined by their ability to take up S. aureus (CTV staining lies within the cell boundary) and induce a ROS response (bright DHR signal) (Rabinovitch, 1995). In zebrafish, professional phagocytes are mainly granulocytes and monocytic cells and can be discriminated from each other based on morphological differences, such as cell size, granularity, and nuclear shape (Wittamer et al., 2011). To compare samples incubated with CTV-S. aureus and the samples where phagocytosis is inhibited (CTV-S. aureus + CCB and CTV-S. aureus + Ice), we used the statistical analysis included in Image3C based on a negative binomial regression model (Figure 3B, C, Figure 3—figure supplement 2; Supplementary files 5 and 6). Statistical analyses reported clusters with differences in relative abundance between phagocytosis and phagocytosis-inhibited samples. Visualizing these clustered event images (Supplementary file 7) while considering the values and intensities of their morphological and fluorescent features (Supplementary file 3) allowed identification of 4 clusters of professional phagocytes: granulocytes within clusters Dr4_P, Dr12_P, and Dr13_P and monocytic cells in cluster Dr21_P (Figure 3A, B). The morphology of cells in cluster Dr12_P is characteristic of phagocytic neutrophils (Figure 2B, Figure 3A) that become adhesive and produce extracellular traps upon recognition of bacterial antigens (Palić et al., 2007). Overall, the relative abundance of professional phagocytes is 5–10% (Figure 3C), which is in line with previous studies that estimated the number of professional phagocytes in WKM tissue of adult zebrafish using classical morphological approaches (Wittamer et al., 2011). It is also noteworthy that in line with other studies (Page et al., 2013) we did not observe a cluster of lymphocytes (e.g., B-cells) that actively phagocytize CTV-S. aureus bacteria (Figure 2; Supplementary file 7). Compared to the classical morphological approaches, Image3C allows to analyze thousands of events in a high-throughput and unbiased fashion, allowing the study of rare cell morphologies and increasing results confidence and reproducibility. These results show that Image3C can successfully analyze biological processes since we were able to recapitulate the presence, cell type, and frequency of professional phagocytes in adult zebrafish WKM organ.

Figure 3 with 2 supplements see all

Download asset Open asset

A new aspect that Image3C highlighted is that CCB selectively affects cell viability based on cell identity, introducing artifacts and cell damage, actions not specific to inhibition of phagocytosis (Figure 3B). All mature erythrocyte-containing clusters had a significantly higher cell count in the CTV-S. aureus samples compared to the CTV-S. aureus + CCB ones (Figure 3B; Supplementary files 3 and 5). Cluster analysis revealed that erythrocytes were almost absent in samples incubated with CCB (Supplementary file 3), while there was a significant increase in the relative abundance of clusters containing dead and apoptotic cells (Figure 3B; Supplementary file 5). Both outcomes are likely due to reduced cell viability of erythrocytes upon CCB incubation. Moreover, we excluded the possibility of higher cell death in the professional phagocytes upon CCB incubation since pseudo-phagocytes (phagocytes with DHR response but no internalized CTV-S. aureus) were significantly more abundant in the CTV-S. aureus + CCB sample (Figure 3B; Supplementary file 5). These results are remarkable since Image3C allowed us to observe a specific effect of CCB on erythrocytes' viability in zebrafish that, as far as we know, was not described before.

Image3C analysis also uncovered another important biological observation. When we inhibited phagocytosis by incubating the single cell suspension on ice (CTV-S. aureus + Ice) and compared the specificity of inhibition with the CTV-S. aureus + CCB sample (Figure 3C; Supplementary file 6), we discovered that the inhibition of phagocytosis through low temperature only affects adhesive neutrophils (cluster Dr12_P) (Figure 3C). This is suspected to occur as ice inhibits adhesion, while CCB effectively blocks phagocytosis in all professional phagocytes in zebrafish WKM tissue by acting on the cytoskeleton. The use of Image3C allowed us to specifically identify cell types that are sensible to low temperature and those that are not, confirming the existence of different phagocytosis mechanisms and providing additional knowledge about pros and cons of different protocols that can be applied to inhibit phagocytosis based on specific goals and needs.

Since we aimed to provide a tool that is widely applicable, we tested Image3C versatility on the apple snail Pomacea canaliculata, an emerging organism for which molecular and cell biological tools have yet to be fully developed. As such, we repeated the same experiments done in zebrafish on the hemolymph of P. canaliculata. For morphological examination of the cellular composition of the hemolymph collected from female adults in homeostasis conditions, we stained the single cell suspensions with Draq5 (DNA dye) and ran on the ImageStream^X Mark II. We used Image3C to analyze the images of the events and identified 9 cell clusters (Figure 4A, Figure 4—figure supplement 1). Two of these clusters comprised cell doublets, debris, and dead cells (clusters Pc5 and Pc8) and the other clusters, based on inspection of cell images, were grouped into two main categories (Figure 4A; Supplementary file 8). The first category includes small blast-like cells (cluster Pc4) and intermediate cells (clusters Pc2 and Pc3) with high nuclear-cytoplasmic (N/C) ratio. These cells morphologically resemble the Group I hemocytes previously described using a classical morphological approach (Accorsi et al., 2013). The second category comprised larger cells with lower N/C ratio and abundant membrane protrusions (clusters Pc1, Pc6, Pc7, and Pc9). Likely, these cells correspond to the previously described Group II hemocytes that include both granular and agranular cells (Accorsi et al., 2013).

Figure 4 with 1 supplement see all

Download asset Open asset

To identify which of these clusters were enriched for granular cells, we looked at the heatmap with feature values for each individual cluster (Figure 4B; see Supplementary file 1 for feature description). Cluster Pc6 had the highest values for the features related to cytoplasm texture and granularity (i.e., area granularity, intensity granularity, and signal granularity) amongst all clusters other than cell doublets (Figure 4B; Supplementary files 3 and 8). Based on these data, we identified cluster Pc6 as the one containing the granular hemocytes. The clusters obtained by Image3C were not only homogeneous and biologically meaningful but were also consistent with published P. canaliculata hemocyte classification obtained by classical morphological methods (Accorsi et al., 2013). Such remarkable consistency was observed both in terms of identified cell morphologies and their relative abundance in the population of circulating hemocytes (Figure 4C; Supplementary file 8). For example, the relative abundance of the previously reported small blast-like cells is 14.0%, a value almost identical to the abundance of the corresponding cluster Pc4 (13.8%).

Similarly, the category of larger hemocytes or Group II hemocytes represents 80.4% of the circulating cells as measured by traditional morphological methods, while clusters Pc1, Pc6, Pc7, and Pc9 combined represent 72.4% of the events analyzed with Image3C (Figure 4C; Supplementary file 3). A subset of these cells are the granular cells (cluster Pc6), which correspond to 7.7% of all hemocytes by classical histological methods and 8.9% by Image3C. The intermediate cells (clusters Pc2 and Pc3) are the least represented in both approaches, with a relative abundance of 5.6% and 10.6% for the manually and Image3C analyzed events, respectively (Figure 4C; Supplementary file 3). This difference is likely best explained by the remarkable difference in both the number of cells and the number of features that can be considered for analysis by Image3C. Only a few hundred hemocytes were visually analyzed using traditional histological methods based only on cell diameter and N/C ratio (Accorsi et al., 2013). In contrast, the automated pipeline used by Image3C facilitated the analysis of 10,000 nucleated events for each sample and considered 25 morphological features for each cell. The significantly higher number of morphological features simultaneously considered by Image3C also explains the higher number of clusters and improved resolution to distinguish cell types compared to the traditional methods. Hence, Image3C not only can properly analyze cells obtained from an emerging research organism generating biologically meaningful and informative clusters but also represents an unprecedented increase in the accuracy of cell type identification over traditional histological methods, while also allowing high-throughput capability.

As with zebrafish, we also performed a phagocytosis experiment on hemocytes from P. canaliculata. Our goal was to test if it is possible with an emerging research organism to successfully discover cell phenotypes and functions and obtain information about specific biological processes of interest by using Image3C to compare cell populations among treated and control samples.

Here, we setup the phagocytosis assay incubating the cells with CTV-S. aureus and DHR at room temperature. The phagocytosis was inhibited, as control, either adding EDTA (CTV-S. aureus + EDTA) or using low temperature by incubating samples on ice (CTV-S. aureus + Ice). Events collected on the ImageStream^X Mark II were analyzed with Image3C and clustered in 20 distinct clusters using quantifications of morphological and fluorescent features (see Supplementary file 2 for feature description), including nuclear staining, phagocytized S. aureus, and DHR positivity (Figure 5A, Figure 5—figure supplement 1). We compared the phagocytosis-permissive samples (CTV-S. aureus) with samples where phagocytosis was inhibited by EDTA incubation or low temperature using the statistical analysis included in Image3C based on a negative binomial regression model (Figure 5B, C, Figure 5—figure supplement 2; Supplementary files 9 and 10). The clusters with relative abundance significantly higher in the phagocytosis samples (Figure 5B; Supplementary files 3 and 11) and with high intensities of both DHR and bacteria signals (Figure 5—figure supplements 3 and 4) are the two clusters that we identify as enriched with professional phagocyte (cluster Pc5_P and Pc17_P) (Figure 5B, Figure 4—figure supplement 1, Figure 5—figure supplement 4; Supplementary file 11). The two clusters show a different DHR signal intensity (ROS response) from one another upon bacteria exposure (cluster Pc5_P with high DHR signal, cluster Pc17_P with low DHR signal) (Figure 5—figure supplement 3; Supplementary files 3 and 11). Both Pc5_P and Pc17_P relative abundance is significantly higher in the phagocytosis samples compared to the EDTA-treated sample (Figure 5C; Supplementary file 9), showing that EDTA effectively inhibits phagocytosis for both types of professional phagocytes. In the sample where the phagocytosis was inhibited by low temperature, however, only cluster Pc17_P had a significantly lower relative abundance compared to the phagocytosis sample (Figure 5C, Figure 3—figure supplement 2; Supplementary file 10). We can conclude that similar to CCB inhibition in the zebrafish phagocytosis experiment, EDTA is a more effective and generalized (not cell type-specific) inhibitor of phagocytosis than low temperature. These results show that also in an emerging research organism Image3C allowed discovery of new aspects of this biological process and highlighted differences among professional phagocytes that would have been difficult to detect with other methods.

Figure 5 with 4 supplements see all

Download asset Open asset

The data analysis with Image3C clearly highlighted that CCB and EDTA, two classical phagocytic inhibitors commonly used in controls for phagocytosis experiments in vertebrates and invertebrates, respectively, result in a drastic change of cell morphology and cell viability. This consequence is not easily detectable by other methods and is therefore often overlooked. In the present work, these changes significantly modified the overall cell cluster number and distribution and indicate that the effects of CCB and EDTA on cell morphology should be taken into consideration in any study of morphological features of cells with phagocytosis properties because artifacts might be significant.

When determining differences between control and experimental treatments, Image3C necessarily combines images and data from all samples and then clusters the cells. This must be taken into consideration for experimental planning. Experiments meant to analyze cell composition and morphological diversity in one biological domain (e.g., homeostasis condition) (Figure 2, Figure 4) should be carried out separately from those in other domains that are likely to introduce changes in the cell population composition or cell morphologies representing a confounding factor for the de novo clustering in homeostasis condition. Image3C clustering works best when used, at the same time, only on samples belonging to a single experimental domain, such as homeostasis or the phagocytosis assay. An issue that emerges when analyzing different experimental sets independently is the increase of time requirement for analytical steps, the likelihood of introducing errors, and the need to repeatedly annotate the clusters in the FDL graph obtained from each experimental set. This last element is required for comparing cell type behaviors among multiple experiments and provide a global understanding of their functions and response to treatments (i.e., cluster #1 from one analytical run cannot be expected to match cell morphologies with cluster #1 from another run since there is a stochastic element to the process).

This last point is probably the most challenging since mistakes can easily be introduced based on user biases or lack of sufficient pre-existing knowledge about cell morphologies or of cell biomarkers that would allow a confident cross-annotation between multiple FDL graphs. In addition, we observed that the number of clusters drastically increases when including treatments that influence cell morphological properties of the cell. As an example, while we detected 9 unique clusters in naïve hemolymph samples, we detected 20 clusters in the phagocytosis experiment (Figure 3A, Figure 4A). This is in part due to the fact that professional phagocytes change their morphology upon detection of pathogens (Palić et al., 2007), thus creating new clusters. Similarly, the complexity of the clustering is also increased by treatments, such as CCB and EDTA incubations, that are necessary to ensure identification of professional phagocytes, but have a strong impact on the morphology of the cells making the clustering and annotation steps more challenging and prone to mistakes since treated samples contain aberrant populations not found at homeostasis (Figure 5A; Supplementary file 11).

To provide an alternative for streamlining the analysis of multiple experimental sets upon initial de novo clustering and cell type identification in homeostasis samples, we included in Image3C the possibility to use these initial images and their cluster IDs to train a CNN without manually classifying the images (Figure 1). This trained classifier can then be used to assign the cell images subsequently collected from additional experimental sets to one of the clusters defined in the homeostasis condition in a high-throughput way. In this way, it will be possible to determine the behavior of a specific cell type through multiple experimental sets without re-clustering whenever new data is acquired. A crucial element that allows this approach is also represented by the ImageStream^X Mark II system that provides highly reproducible and comparable images of cells coming from different experiments and acquired at different days, introducing much less variability than standard light or electron microscopy.

For our pipeline, then, a CNN (LeCun et al., 1989) based on the architecture of DenseNet (Huang et al., 2017) was deployed to (1) use, as training set, images and clusters obtained from a first group of samples (e.g., homeostasis conditions, naïve cells, or WT samples) analyzed in an unbiased way by de novo clustering and (2) assign new cell images acquired through ImageStream^X Mark II system to their respective classes. As proof of concept, we used the clusters identified for P. canaliculata hemocytes in homeostasis condition with the first part of the pipeline (Figure 4A) for training and setting up the CNN classifier. This approach would define the classes based on the unbiased de novo clustering of thousands of cells with no need for formal annotation or previous knowledge about cell types and tissue composition. To prepare a dataset for training the classifier, we first combined clusters that strongly overlapped with one another in terms of morphological characteristics (e.g., doublets and dead cells) to increase accuracy of the classifier (Figure 6A). We used 80% of the cells obtained in the original P. canaliculata dataset together with their cluster IDs to train the classifier through over 25,000 iterations. After each iteration, we tested the training with 10% of the original dataset and determined the relative accuracy by scoring numbers of cells whose cluster ID assigned by the classifier matched the original cluster ID (Figure 6B, C). The remaining 10% of the original dataset was used to calculate the precision of the trained classifier. Clusters with higher support numbers obtained higher precision scores. The weighted average precision score (f1-score, precision average score across clusters controlling for support numbers) of 0.75 is relatively high considering the complexity of the phenotype (brightfield [BF], darkfield, and Draq5 images) (Figure 6D) and comparable to other studies using ML for cell classification (Blasi et al., 2016). The true probability match for each individual cell (probability for each cell that the class assigned by the classifier would match the original cluster ID) demonstrated that lower true probability matches occurred where clusters strongly overlapped or where cell phenotypes are intermediate between clusters, providing an additional layer of information about our dataset (Figure 6D).

Figure 6

Download asset Open asset

To test the efficiency of this pipeline, we extracted all the images belonging to the two clusters identified in the phagocytosis assay as cluster-containing phagocytes and determined to which naïve cell type they correspond using the CNN classifier and only the BF, SSC, and Draq5 channels (i.e., DHR and labeled bacteria signals were not used). We found that 59.4, 6.2, and 9.2% of the phagocytes belonged to cluster Pc1_CT, Pc6_CT, and Pc7_CT, respectively (Figure 6E), where CT stands for classifier training. These results confirmed a previously published result that used classical morphological staining and manual annotation to conclude that the hemocytes able to phagocytize were primarily Group II hemocytes (Accorsi et al., 2013). Only 8% of the phagocytes were clustered in the Group I hemocytes, here represented by clusters Pc2_CT, Pc3_CT, and Pc4_CT, while the remaining 17.2% were assigned by the CNN to the cluster Pc5_CT, constituted by doublets and dead cells (Figure 6E). This result can be explained by the fact that in vitro phagocytosis triggers microaggregate formation (hemocyte-hemocyte adhesion) in invertebrate hemocytes that resembles the nodule formation observed in vivo (Walters, 1970). It is important to observe how this analysis allowed us to assign phagocytes to cell types using the annotation already performed in Figure 4A (de novo clustering of hemocytes in homeostasis condition) without the need to reannotate the FDL obtained during the phagocytosis assay (Figure 5A).

To test the adaptability of the trained CNN to new datasets, we collected hemocytes from male apple snail specimens, stained the cells with Draq5, and recorded BF, SSC, and nuclei images from 10,000 cells on the ImageStream^X Mark II as previously described. We extracted the images of the cells and used the trained CNN classifier to determine the relative abundance of hemocytes collected from male snails in the seven classes of the classifier (Figure 6F). First, we visually compared the female hemocytes clustered by the de novo clustering with the male hemocytes that were run on the ImageStream^X Mark II and were assigned to a class by the classifier (Figure 6F). This comparison shows that the female and male hemocytes belonging to the same cluster are morphologically extremely similar and different from the hemocytes assigned to other clusters (Figure 6F). This demonstrates that the CNN classifier can be trained with a first group of samples and then it can successfully analyze new datasets acquired later on. The comparison between female and male hemocyte compositions revealed that the clusters significantly different in terms of relative abundance are Pc1_CT and Pc7_CT (Group II agranular large hemocytes) and Pc6_CT (Group II granular large hemocytes) (Figure 6F). Significantly, prior studies detected no differences between females and males hemocytes composition through manual classification and counting using a classical morphological approach (Accorsi et al., 2013). The reduced user bias and high-throughput analysis presented here, in contrast, allowed us to determine that one of the two subpopulations of agranular large hemocytes was significantly more abundant in the female animals (Pc1_CT: 53 and 38% in females and males, respectively) while the other agranular (Pc7_CT) as well as the granular large hemocytes (Pc6_CT) was significantly more abundant in the male animals (Pc7_CT: 14 and 20% in females and males, respectively; Pc6_CT: 6 and 10% in females and males, respectively) (Figure 6F).

While the biological significance of this observation is not going to be further investigated in this paper, the discovery highlights the power of Image3C analysis compared to traditional methods for determining and quantifying the composition of cell populations. These experiments demonstrate that Image3C, in combination with the presented convolutional classifier, can analyze large experimental datasets and identify significances with small effect sizes. Importantly, Image3C analysis is independent of observer biases and does not require prior knowledge about expected tissue composition or the expected effect of treatment on cell morphology.

We have developed a powerful new method to analyze at single-cell resolution the composition of any cell population obtained from research organisms for which species-specific reagents (such as fluorescently tagged antibodies), biomarkers, single-cell atlases, or a high-quality genome for a scRNA-seq approach are not available. We demonstrated that Image3C can identify different cell populations based on morphology and/or function through de novo clustering and highlight important changes in cell type abundance and cell population composition caused by experimental or natural perturbation (sex, treatment, experimental protocol). Image3C does not require, at any point, prior knowledge about the tissue composition or cell type-specific markers, although, if available, they can be included and used. Furthermore, in combination with the CNN classifier trained on these clusters, we demonstrate that Image3C is capable of bias-free and high-throughput analysis of large experimental datasets making it possible to compare a specific cell type behavior or population composition across multiple experiments. Image3C is extremely versatile and can be applied to any tissue or cell population of interest and is adaptable to a variety of experimental designs. Although Image3C was developed in response to the need of analyzing cell composition of tissues in emerging research organisms, the Image3C tool could be potentially used also to add to transcriptomic dataset an additional and complementary layer of information based on cell morphology. Given the recent advancement in image-based flow cytometry that enables image capturing together with cell sorting (Nitta et al., 2018), a scRNA-seq approach in combination with the Image3C pipeline would enable simultaneous analysis of both the morphological/phenotypic and genetic properties of a cell population at single-cell resolution.

Twelve-month-old, wild type, female, adult zebrafish were euthanized with cold 500 mg/L MS-222 solution for 5 min. WKM was dissected as previously described (Traver et al., 2003) and transferred to 40 µm cell strainer with 1 mL of L-15 media containing 10% water, 10 mM HEPES, and 20 U/mL Heparin (L-90). Cells were gently forced through the cell strainer with the plunger of a 3 mL disposable syringe. The strainer was washed once with 1 mL of L-90 and the resulting single cell suspension was centrifuged at 500 rcf at 4 °C for 5 min. The supernatant was discarded, and the cells were resuspended in 1 mL of L-90 containing 5% fetal calf serum (FCS), 4 mM L-glutamine, and 10,000 U of both penicillin and streptomycin (L-90 media). The cells were counted in a 1:20 dilution on the EC-800 flow cytometer (Sony) using scatter properties.

Specimens of the apple snail P. canaliculata (Mollusca, Gastropoda, Ampullariidae) were maintained and bred in captivity, in a water recirculation system filled with artificial freshwater (2.7 mM CaCl₂, 0.8 mM MgSO₄, 1.8 mM NaHCO₃, 1:5000 Remineralize Balanced Minerals in Liquid Form [Brightwell Aquatics]). The snails were fed twice a week and kept in a 10:14 light:dark cycle. Wild-type adult snails, 7–9 months old and with a shell size of 45–60 mm, were starved for 5 days before the hemolymph collection (Accorsi et al., 2013). If not differently specified, female snails were used for the experiments. The withdrawal was performed by applying a pressure on the operculum and dropping the hemolymph directly into an ice-cold tube. The hemolymph collected from different animals was not pooled together. The hemolymph was immediately diluted 1:4 in Bge medium + 10% fetal bovine serum (FBS) and then centrifuged at 500 rcf for 5 min. The pellet of cells was resuspended in 100 µL of Bge medium + 10% FBS. The Bge medium (also known as Biomphalaria glabrata embryonic cell line medium) is constituted by 22% (v/v) Schneider’s Drosophila Medium, 4.5 g/L lactalbumin hydrolysate, 1.3 g/L galactose, 0.02 g/L gentamycin in MilliQ water, pH 7.0.

WKM cells from zebrafish were isolated as described before and plated at 4 × 10⁵ cells/well in a 96-well plate in 200 µL of L-90 media and incubated for 3 h at room temperature. Cells were stained with 5 µM Draq5 (Thermo Fisher Scientific) for 10 min and subsequently run on the ImageStream^X Mark II (Amnis Millipore Sigma), where 10,000 nucleated and focused events were recorded for each sample (n = 8). Erythrocytes were outgated to enrich for immune-relevant cells and prevent overclustering in the subsequent analysis. The latter is due to the fact that fish erythrocytes are nucleated and their biconcave shape results in different morphological feature values only depending on their orientation during image acquisition.

The P. canaliculata hemocytes were stained with 5 µM Draq5 for 10 min, moved to ice, and subsequently run on the ImageStream^X Mark II, where 10,000 nucleated and focused events were imaged for each sample (n = 5).

Staphylococcus aureus (Thermo Fisher Scientific) were resuspended in PBS at the final concentration of 100 mg/mL and incubated with 5 µM CTV (Thermo Fisher Scientific) for 20 min. Labeled bacteria were centrifuged and resuspended in PBS for three times to remove unbound dye and then stored at −20 °C as single-use aliquots. Cells, obtained from fish WKM or snail hemolymph and in a single cell suspension, were plated in a 96-well plate at a concentration of 4 × 10⁵ cells/well in 200 µL of medium and incubated with 2 × 10⁷ CTV-coupled S. aureus/well for 3 h at room temperature. As control, the phagocytosis was inhibited incubating the cells + CTV-S. aureus mix either on ice (for both species) or with 0.08 mg/mL CCB for zebrafish cells or with 30 mM EDTA and 10 mM HEPES for apple snail cells (Cueto et al., 2015; Li et al., 2006). After 2 h and 30 min, we added 5 µM DHR (Thermo Fisher Scientific) to the cell suspension to stain cells positive for ROS production. To control for this treatment with DHR, we incubated one aliquot of cells with 10 ng/mL phorbol 12-myristate 13-acetate (PMA) to induce ROS production. At 2 h and 50 min since the beginning of incubation with CTV-S. aureus, all the samples were stained with 5 µM Draq5 for 10 min. After 3 h incubation with bacteria, cells were moved and stored on ice and subsequently run on the ImageStream^X Mark II, where 10,000 nucleated and focused events were imaged for each sample (at least n = 4 snail and n = 6 fish) at a speed of 1,000 images/s.

Following cell preparation, data were acquired from each sample on the ImageStream^X Mark II at 60× magnification, slow/sensitive flow speed (1,000 images/s), using 633, 488, and 405 nm laser excitation. BF was acquired on channels 1 and 9, DHR (488 nm excitation) on channel 2, CTV-S. aureus (405 nm excitation) on channel 7, Draq5 (633 nm excitation) on channel 11, and SSC was acquired on channel 6. Single-color controls were also acquired for each fluorescent channel to allow for fluorescence spillover correction.

An interactive map representing the pipeline, the software used, the format of the exported files, and an approximation of time required for running the individual steps is provided in Figure 1—figure supplement 1. Raw images from the ImageStream^X Mark II system (RIF files, a type of modified 16-bit TIFF file) were compensated (spillover and other corrections applied), background was subtracted, and features were calculated using IDEAS 6.2 software (Amnis Millipore, free for download once an Amnis user account is created). The resulting compensated image files (CIF files) were used to quantify features for all cells and samples. Supplementary files 1 and 2 report the list of features used for each organism and for each experiment and their description. These per-object feature matrices (DAF files) were then exported from IDEAS into FCS files. Exported FCS files were processed in R (R Development Core Team, 2014). In order to trim redundant features that contribute noise but little new information, Spearman’s correlation values for each pair of features were calculated using all events of a representative sample and one of the features of the pair was trimmed when correlation between the two was ≥0.85 (Figure 2—figure supplement 1; Caicedo et al., 2017). The Spearman’s correlation of the mean values of remaining features per each sample was then used to identify outliers among sample replicates. Samples with correlation of mean feature values below 0.85 with the set were discarded (Figure 2—figure supplement 2, Figure 3—figure supplement 1, Figure 4—figure supplement 1, Figure 5—figure supplement 1), although in general the replicates were consistent. Also, while morphological features did not require any transformation or normalization, fluorescence intensity features were transformed using the estimateLogicle() and transform() functions from the R flowCore package (Ellis et al., 2018; Hahne et al., 2009) to improve homoscedasticity (homogeneity of variance) of distributions. DNA intensity features were also normalized to align all 2N and 4N peak positions and remove intensity drift between samples (Figure 2—figure supplement 3) using the gaussNorm() function from flowStats package (Hahne et al., 2018). The processed data was exported from R (R Development Core Team, 2014) using writeflowSet() function in flowCore package (Ellis et al., 2018; Hahne et al., 2009) as CSV or FCS files, depending on downstream needs for the file output.

These processed data files were then imported into the VorteX clustering environment for X-Shift k-nearest-neighbor clustering (free to install) (Samusik et al., 2016). X-Shift was selected as a clustering method for Image3C based on a previously published analysis and comparison of clustering methods (Weber and Robinson, 2016). From that work, we determined that X-shift represents an optimal trade-off: identifying low-frequency populations, accurately identifying ‘true’ clusters (i.e., F1 scores), not requiring for a priori knowledge of the number of clusters (populations) and having reasonable runtimes (due to hardware CPU requirements). We also did an early comparison of X-shift with K-means (data not shown) and determined that K-means was insufficient for our purposes as known cell populations were not well represented by clusters, and we did not want to specify the number of expected clusters into the method’s input parameters since this would not be known in experimental use. During the import into VorteX, all features were scaled to 1 SD to equalize the contribution of features towards clustering. Clustering was performed in VorteX testing a range of k values (typically from 5 to 150), choosing a final k value using the ‘find elbow point for cluster number’ function in VorteX and confirming visually that over- or underclustering did not occur. FDL graphs of a subset of cells obtained from each set of samples were also generated in VorteX, and cell coordinates in the resultant 2D space were exported along with graphML representation of the FDL graph. Finally, tabular data (CSV files) was exported from VorteX including a master table of every cell event with its cluster assignment and original sample ID, as well as a table of the average feature values for each cluster and counts of cells per cluster and per sample.

Clustering results were further analyzed and plotted in R (R Development Core Team, 2014) by merging all cell events and feature values with cluster assignments and X/Y coordinates for FDL graph. Using this merged data and the graphML file exported from VorteX, new FDL graphs were created for each treatment condition using the igraph package (Csardi and Nepusz, 2006) in R (R Development Core Team, 2014). Statistical analysis of differences in cell counts per cluster by condition was performed using negative binomial regression of cell counts per cluster, plots of statistic results and other results were generated, and CSV files containing cell ID, sample ID, feature values, and X/Y coordinates in FDL graph were exported for each sample. The subsequent use of FCS Express Plus version 6 (DeNovo software, free alternative are mentioned later in the text) allowed visualization of cell images using DAF/CIF files by cluster and customized subsets of the FDL graphs.

DAF files were opened in FCS Express Plus software, and the ‘R add parameters’ transformation feature with a custom script was used to merge the clustering data saved in the CSV files generated above with both DAF and CIF files (feature values and image sets, respectively). FCS Express Plus was utilized at this stage of work because it is the only platform currently available that works with Amnis DAF and CIF files while also running transformation processes driven by R scripting. ImageJ Bio-Formats allows reading images from DAF and CIF files, but we got pixels with a value much higher than expected, probably due to a bug that has not been fixed yet. This allowed to visualize image galleries of cells within each cluster and gate by features of interest on 2D plots (more traditional flow cytometry analysis) for exploring the clustering results and identifying clusters and populations of interest. FCS Express Plus is a proprietary software, but similar results can be obtained with IDEAS software where a text file with cluster IDs for each image event can be imported and the cluster information can be matched to the event images.

The full complement of R packages used includes flowCore (Ellis et al., 2018; Hahne et al., 2009), flowStats (Hahne et al., 2018), igraph (Csardi and Nepusz, 2006), ggcyto (Jiang, 2015), ggridges (Wilke, 2018), ggplot2 (Wickham, 2016), stringr (Wickham, 2010), hmisc (Harrell and Dupont, 2019), and caret (Kuhn, 2008). Figure 1—figure supplement 1 can be used as an interactive map of Image3C pipeline, where, upon clicking on the different portions of the pipeline, the users will be automatically directed to the corresponding sections of our GitHub repository. The GitHub repository at https://github.com/stowersinstitute/LIBPB-1390-Image3C reports a detailed description of all these processing steps and includes sample scripts, workflow files, and example datasets and tutorials.

Once the clusters were defined with the previously described de novo clustering analysis, we used a CNN (LeCun et al., 1989) based on the architecture of DenseNet (Huang et al., 2017) for image classification. Out of all the events captured with the ImageStream^X Mark II system, we selected only single nucleated objects applying gates on area vs. aspect ratio plot and Draq5 intensity plot to achieve this selection, respectively. For these objects, we exported 16-bit TIFF images (one channel per fluorescence/BF image ‘color’) using IDEAS 6.2 software.

Because images from the ImageStream^X Mark II system have nonuniform sizes, each image was cropped or padded to 32 × 32 pixels using NumPy indexing (van der Walt et al., 2011) in a Python script. The CNN consists of three dense blocks that transition from three-channel image input of 32 × 32 × 3 to a final size of 4 × 4 × 87 with 87 feature maps. A dense block includes three convolution layers, each followed by leaky ReLU activation. The last step of the block is a strided convolution used to down-sample the width and height of the feature maps by a factor of 2. The final layer of the CNN flattens the 4 × 4 × 87 array into a 1D vector of length 1392 and is fully connected to the output layer, that is, a 1D vector with a length of the number of classes for prediction. The CNN used softmax cross-entropy for the loss function with L2 regularization, and the Adam optimizer (Kingma and Ba, 2014) was used to minimize the loss function. The CNN was implemented in Python using the TensorFlow platform (Abadi et al., 2015) and the SciPy ecosystem (Oliphant, 2006; Oliphant, 2007; Pedregosa et al., 2011).

The CNN was used to train a classifier using over 35,000 images of P. canaliculata hemolymph cell types in homeostasis condition acquired with the ImageStream^X Mark II system. The event images were randomly split in 80% for training, 10% for testing during training, and 10% for final validation. The first 80% of the images were used together with their cluster IDs obtained by the de novo clustering to train the classifier. The learning rate for the Adam optimizer was set to 0.0006 with a decaying learning rate starting at 0.001 and decreasing by 1% each step. The training proceeded for 25,000 iterations with a size of 256 randomly selected images for each iteration. After each iteration, 10% of the cells of the original P. canaliculata dataset was used to test the classifier. The loss and accuracy of the CNN were recorded after every 100 iterations to monitor the performance. The CNN loss was defined by the softmax of the cross-entropy (Dahal, 2017) between the final output and the one-hot-encoded image labels. To avoid the CNN memorizing the training set, L2 regularization was applied to the weights. The training and test sets follow the same accuracy and loss trends over all iterations, indicating the training set is not memorized and can generalize to predict the test set.

The finally trained classifier was tested on the remaining 10% of images that were completely new for the CNN. The trained model was saved for future use, so new images can be inferred by the network to predict cell types. The inference is very fast because only one forward pass is made through the network and no backpropagation occurs. The result of the inference is a vector with length equal to the number of cell type classes. Each element of the vector will be the probability of the cell belonging to the corresponding class; the sum of the vector must be 1. Inferring a complete experiment will provide a probability vector for each image; the list of probability vectors can be saved as CSV file. For new images, the inference results will need to be examined to ensure the predictions are reliable. A large majority of the probability vectors should have a maximum greater than 0.5, and a subset of the images should be visually inspected to verify proper class assignment.

The interactive map of Image3C pipeline (Figure 1—figure supplement 1) includes also the training and the use of the CNN. The GitHub repository at https://github.com/stowersinstitute/LIBPB-1390-Image3C reports a detailed description of all these processing steps, includes sample scripts, workflow files, and example datasets and tutorials and can be easily accessed by clicking on the right side of the Image3C interactive map (Figure 1—figure supplement 1).

Negative binomial regression was performed on tables of cell counts per cluster and per sample and plots were generated using R (R Development Core Team, 2014) with the edgeR package (Robinson et al., 2010), which was developed for RNA-seq analysis, but includes generally applicable and user-friendly wrappers for regression, modeling analysis, and plotting of results. For the comparison of cellular hemolymph composition between females and males of P. canaliculata, a one-way ANOVA with subsequent FDR (Benjamini–Hochberg, correction for multiple testing) was used.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

We acknowledge Hua Li for her assistance on the statistical analysis and also thank the Laboratory Animal Services and the Aquatics Facility at the Stowers Institute for Medical Research for animal husbandry. We would like to thank Blair Benham-Pyle, Carolyn Brewster, Viraj Doddihal, Julia Peloggia de Castro, and Barbara Milutinović for their critical comments on an earlier version of this manuscript. This work was supported by institutional funding to ACB, CW, ASA, and NR. ASA is a Howard Hughes Medical Institute Investigator. NR is further supported by the Edward Mallinckrodt Foundation, NIH Grant R01 GM127872, DP2AG071466, and NSF EDGE award 1923372. AA was supported by the Emerging Models grant from the Society for Developmental Biology (SDB) and the postdoctoral fellowship from the American Association of Anatomists (AAA). RP was supported by a grant from the Deutsche Forschungsgemeinschaft (PE 2807/1-1).

Animal experimentation: Research and animal care were approved by the Institutional Animal Care and Use Committee (IACUC) of the Stowers Institute for Medical Research. protocol (#2018-0187 and #2019-080).

© 2021, Accorsi et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.