Highly multiplexed immunofluorescence imaging of human tissues and tumors using t-CyCIF and conventional optical microscopes

The architecture of normal and diseased tissues strongly influences the development and progression of disease as well as responsiveness and resistance to therapy. We describe a tissue-based cyclic immunofluorescence (t-CyCIF) method for highly multiplexed immuno-fluorescence imaging of formalin-fixed, paraffin-embedded (FFPE) specimens mounted on glass slides, the most widely used specimens for histopathological diagnosis of cancer and other diseases. t-CyCIF generates up to 60-plex images using an iterative process (a cycle) in which conventional low-plex fluorescence images are repeatedly collected from the same sample and then assembled into a high-dimensional representation. t-CyCIF requires no specialized instruments or reagents and is compatible with super-resolution imaging; we demonstrate its application to quantifying signal transduction cascades, tumor antigens and immune markers in diverse tissues and tumors. The simplicity and adaptability of t-CyCIF makes it an effective method for pre-clinical and clinical research and a natural complement to single-cell genomics.


Introduction
Histopathology is among the most important and widely used methods for diagnosing human disease and studying the development of multicellular organisms. As commonly performed, imaging of formalin-fixed, paraffin-embedded (FFPE) tissue has relatively low dimensionality, primarily comprising Hematoxylin and Eosin (H&E) staining supplemented by immunohistochemistry (IHC). The potential of IHC to aid in diagnosis and prioritization of therapy is well established (Bodenmiller, 2016), but IHC is primarily a single-channel method: imaging multiple antigens usually involves the analysis of sequential tissue slices or harsh stripping protocols (although limited multiplexing is possible using IHC and bright-field imaging [Stack et al., 2014;Tsujikawa et al., 2017]). Antibody detection via formation of a brown diamino-benzidine (DAB) or similar precipitates are also less quantitative than fluorescence (Rimm, 2006). The limitations of IHC are particularly acute when it is necessary to quantify complex cellular states and multiple cell types, such as tumor infiltrating regulatory and cytotoxic T cells (Postow et al., 2015) in parallel with tissue and pharmaco-dynamic markers.
Advances in DNA and RNA profiling have dramatically improved our understanding of oncogenesis and propelled the development of targeted anticancer drugs (Garraway and Lander, 2013). Sequence data are particularly useful when an oncogenic driver is both a drug target and a biomarker of drug response, such as BRAF V600E in melanoma (Chapman et al., 2011) or BCR-ABL in chronic myelogenous leukemia (Druker and Lydon, 2000). However, in the case of drugs that act through cell non-autonomous mechanisms, such as immune checkpoint inhibitors, tumor-drug interaction must be studied in the context of multicellular environments that include both cancer and non-malignant stromal and infiltrating immune cells. Multiple studies have established that these components of the tumor microenvironment strongly influence the initiation, progression and metastasis of cancer (Hanahan and Weinberg, 2011) and the magnitude of responsiveness or resistance to immunotherapies (Tumeh et al., 2014).
Single-cell transcriptome profiling provides a means to dissect tumor ecosystems at a molecular level and quantify cell types and states (Tirosh et al., 2016). However, single-cell sequencing usually requires disaggregation of tissues, resulting in loss of spatial context (Tirosh et al., 2016;Patel et al., 2014). As a consequence, a variety of multiplexed approaches to analyzing tissues have recently been developed with the goal of simultaneously assaying cell identity, state, and morphology (Giesen et al., 2014;Gerdes et al., 2013;Micheva and Smith, 2007;Remark et al., 2016;Gerner et al., 2012). For example, FISSEQ (Lee et al., 2014) enables genome-scale RNA profiling of tissues at single-cell resolution, and multiplexed ion beam imaging (MIBI) and imaging mass cytometry achieve a high degree of multiplexing using antibodies as reagents, metals as labels and eLife digest To diagnose a disease such as cancer, doctors sometimes take small tissue samples called biopsies from the affected area. These biopsies are then thinly sliced and treated with dyes to identify healthy and cancerous cells. However, clinicians and scientists often need to look into what happens inside individual cells in the tissues so they can understand how cancers arise and progress. This helps them to identify different types of tumor cells and to tailor the best treatment for the patient.
To do so, a number of proteins (the molecules involved in nearly all life's processes) need to be tracked in healthy and diseased cells and tissues. This can be done thanks to a range of methods known as immunofluorescence microscopy, but following different proteins on the same slice of a sample is difficult. However, a new type of immunofluorescence known as t-CyCIF may be a solution.
With this technique, a fluorescent compound is applied that will bind to a specific protein of interest. A microscope can pick up the light from the compound when the sample is imaged, which reveals the protein's location in the cell or tissue. Then, a substance is used that deactivates the fluorescence signal. After this, another compound that binds to a new type of protein is used, and imaged. This cycle is repeated several times to locate different proteins. Lastly, the individual images are processed and stitched together to reveal the cells and their internal structures.
Here, Lin, Izar et al. showed that t-CyCIF could be used to study biopsies and to obtain images that covered a large area of healthy human tissues and tumors. The technique helped to track over 60 different proteins in normal and tumor tissue samples from human patients. Several sets of experiments showed that t-CyCIF could uncover the molecular mechanisms that are disrupted during cancer, but also reveal the complexity of a single tumor. In fact, as shown with biopsies of brain cancer, cancerous cells in a tumor can be strikingly different, even when they are close to each other. Finally, the method helped to pinpoint which types of immune cells are involved in fighting a kidney tumor. Overall, such information cannot be obtained with conventional methods, yet is crucial for diagnosis and treatment.
Most laboratories can readily use t-CyCIF since the technique is open source and requires equipment that is easily accessible. In fact, the technique should soon be used to assess how well certain drugs help the immune system combat cancer. Ultimately, better use of biopsies is key to customizing cancer care. mass spectrometry as a detection modality (Giesen et al., 2014;Angelo et al., 2014). Despite the potential of these new methods, they require specialized instrumentation and consumables, which is one reason that the great majority of basic and clinical studies still rely on H&E and single-channel IHC staining. Moreover, methods that involve laser ablation of samples such as MIBI inherently have a lower resolution than optical imaging.
Thus, there remains a need for highly multiplexed tissue analysis methods that (i) minimize the requirement for specialized instruments and costly, proprietary reagents, (ii) work with conventionally prepared FFPE tissue specimens collected in clinical practice and research settings, (iii) enable imaging of ca. 50 antigens at subcellular resolution across a wide range of cell and tumor types, (iv) Figure 1. Steps in the t-CyCIF process. (A) Schematic of the cyclic process whereby t-CyCIF images are assembled via multiple rounds of four-color imaging. (B) Image of human tonsil prior to pre-staining and then over the course of three rounds of t-CyCIF. The dashed circle highlights a region with auto-fluorescence in both green and red channels (used for Alexa-488 and Alexa-647, respectively) and corresponds to a strong background signal. With subsequent inactivation and staining cycles (three cycles shown here), this background signal becomes progressively less intense; the phenomenon of decreasing background signal and increasing signal-to-noise ratio as cycle number increases was observed in several staining settings (see also Multi-scale imaging of t-CyCIF specimens. (A) Bright-field H&E image of a metastasectomy specimen that includes a large metastatic melanoma lesion and adjacent benign tissue. The H&E staining was performed after the same specimen had undergone t-CyCIF. (B) Representative t-CyCIF staining of the specimen shown in (A) stitched together using the Ashlar software from 165 successive CyteFinder fields using a 20X/0.8NA objective. (C) One field from (B) at the tumor-normal junction demonstrating staining for S100-postive malignant cells, a-SMA positive stroma, T lymphocytes (positive for CD3, CD4 and CD8), and the proliferation marker phospho-RB (pRB). (D) A melanoma tumor imaged on a GE INCell Analyzer 6000 confocal microscope to demonstrate sub-cellular and sub-organelle structures. This specimen was stained with phospho-Tyrosine (pTyr), Lamin A/ C and p-Aurora A/B/C and imaged with a 60X/0.95NA objective. pTyr is localized in membrane in patches associated with receptor-tyrosine kinase, visible here as red punctate structures. Lamin A/C is a nuclear membrane protein that outlines the vicinity of the cell nucleus in this image. Aurora kinases A/B/C coordinate centromere and centrosome function and are visible in this image bound to chromosomes within a nucleus of a mitotic cell in prophase (yellow arrow). (E) Staining of a melanoma sample using the GE OMX Blaze structured illumination microscope with a 60X/1.42NA objective shows heterogeneity of structural proteins of the nucleus, including as Lamin B and Lamin A/C (indicated by yellow arrows) and part of the nuclear pore complex (NUP98) that measures~120 nm in total size and indirectly allows the visualization of nuclear pores (indicated by non-continuous staining of NUP98). (F) Staining of a patient-derived mouse xenograft breast tumor using the OMX Blaze with a 60x/1.42NA objective shows a spindle in a mitotic cell (beta-tubulin in red) as well as vesicles staining positive for VEGFR2 (in cyan) and punctuate expression of the EGFR in the plasma membrane (in green). DOI: https://doi.org/10.7554/eLife.31657.006 The following figure supplements are available for figure 2: collect data with sufficient throughput that large specimens (several square centimeters) can be imaged and analyzed, (v) generate high-resolution data typical of optical microscopy, and (vi) allow investigators to customize the antibody mix to specific questions or tissue types. Among these requirements the last is particularly critical: at the current early stage of development of high dimensional histology, it is essential that individual research groups be able to test the widest possible range of antibodies and antigens in search of those with the greatest scientific and diagnostic value. This paper describes a method for highly multiplexed fluorescence imaging of tissues, tissuebased cyclic immunofluorescence (t-CyCIF), inspired by a cyclic method first described by Gerdes et al. (2013). t-CyCIF also extends a method we previously described for imaging cells grown in culture (Lin et al., 2015). In its current implementation, t-CyCIF assembles up to 60-plex images of FFPE tissue sections via successive rounds of four-channel imaging. t-CyCIF uses widely available reagents, conventional slide scanners and microscopes, manual or automated slide processing and simple protocols. It can, therefore, be implemented in most research or clinical laboratories on existing equipment. Our data suggest that high-dimensional imaging methods using cyclic immunofluorescence have the potential to become a robust and widely-used complement to singlecell genomics, enabling routine analysis of tissue and cancer morphology and phenotypes at singlecell resolution.
Results t-CyCIF enables multiplexed imaging of FFPE tissue and tumor specimens at subcellular resolution Cyclic immunofluorescence (Gerdes et al., 2013) creates highly multiplexed images using an iterative process (a cycle) in which conventional low-plex fluorescence images are repeatedly collected from the same sample and then assembled into a high-dimensional representation. In the implementation described here, samples~5 mm thick are cut from FFPE blocks, the standard in most histopathology services, followed be dewaxing and antigen retrieval either manually or on automated slide strainers in the usual manner (Shi et al., 2011). To reduce auto-fluorescence and non-specific antibody binding, a cycle of 'pre-staining' is performed; this involves incubating the sample with secondary antibodies followed by fluorophore oxidation in a high pH hydrogen peroxide solution in the presence of light ('fluorophore bleaching'). Subsequent t-CyCIF cycles each involve four steps ( Figure 1A): (i) immuno-staining with antibodies against protein antigens (three antigens per cycle in the implementation described here) (ii) staining with a DNA dye (commonly Hoechst 33342) to mark nuclei and facilitate image registration across cycles (iii) four-channel imaging at low-and high-magnification (iv) fluorophore bleaching followed by a wash step and then another round of immunostaining. In t-CyCIF, the signal-to-noise ratio often increases with cycle number due to progressive reductions in background intensity over the course of multiple rounds of fluorophore bleaching. This effect is visible in Figure 1B as the gradual disappearance of an auto-fluorescent feature (denoted by a dotted white oval and quantified in Figure 1-figure supplement 1; see detailed analysis below). When no more t-CyCIF cycles are to be performed, the specimen is stained with H&E to enable conventional histopathology review. Individual image panels are stitched together and registered across cycles followed by image processing and segmentation to identify cells and other structures. t-CyCIF allows for one cycle of indirect immunofluorescence using secondary antibodies. In all other cycles antibodies are directly conjugated to fluorophores, typically Alexa 488, 555 or 647 (for a description of different modes of CyCIF see Lin et al., 2015). As an alternative to chemical coupling we have tested the Zenon antibody labeling method (Tang et al., 2010) from ThermoFisher in which isotype-specific Fab fragments pre-labeled with fluorophores are bound to primary antibodies to create immune complexes; the immune complexes are then incubated with tissue samples  ( Figure 1-figure supplement 2). This method is effective with 30-40% of the primary antibodies that we have tested and potentially represents a simple way to label a wide range of primary antibodies with different fluorophores.
Imaging of t-CyCIF samples can be performed on a variety of fluorescent microscopes each of which represent a different tradeoff between data acquisition time, image resolution and sensitivity ( Table 1). Greater resolution (a higher numerical aperture objective lens) typically corresponds to a smaller field of view and thus, longer acquisition time for large specimens. Imaging of specimens several square centimeters in area at a resolution of~1 mm is routinely performed on microscopes specialized for scanning slides (slide scanners); we use a CyteFinder system from RareCyte (Seattle WA) configured with 10 Â 0.3 NA and 40 Â 0.6 NA objectives but have tested scanners from Leica, Nikon and other manufacturers. Figure 2A-B show an H&E image of a~10 Â 11 mm metastatic melanoma specimen and a t-CyCIF image assembled from 165 individual image tiles. The assembly process involves stitching sequential image tiles from a single t-CyCIF cycle into one large image panel, flat-fielding to correct for uneven illumination and registration of images from successive t-CyCIF cycles to each other; these procedures were performed using ImageJ, ASHLAR, and BaSiC software as described in materials and methods (Peng et al., 2017).
In the t-CyCIF image ( Figure 2B) tumor cells staining positive for S100 (a melanoma marker in green [Henze et al., 1997]) are surrounded by CD45-positive immune cells (CD45RO + cells in white) and by stromal cells expressing the alpha isoform of smooth muscle actin (a-SMA in red). By zooming in on one tile, single cells can be identified and characterized ( Figure 2C); in this image, CD4 + and CD8 + T-lymphocytes and proliferating pRB + positive cells are visible. At 60X resolution on a confocal GE INCell Analyzer 6000, kinetochores stain positive for the phosphorylated form of the Aurora A/B/C kinase and can be counted in a mitotic cell (yellow arrowhead in Figure 2D). Nominally superresolution imaging on a GE OMX Blaze Structured Illumination Microscope (Carlton et al., 2010) (using a 60 Â 1.42 Plan Apo objective) reveals very fine structural details including differential expression of Lamin isotypes (in a melanoma, Figure 2E and These data show that t-CyCIF images have readily interpretable features at the scale of an entire tumor, individual tumor cells and subcellular structures. Little subcellular (or super-resolution) imaging of clinical FFPE specimens has been reported to date (but see Chen et al., 2015), but fine subcellular morphology has the potential to provide dramatically greater information than simple integration of antibody intensities across whole cells.
To date, we have tested commercial antibodies against~200 different proteins for their compatibility with t-CyCIF; these include lineage makers, cytoskeletal proteins, cell cycle regulators, the phosphorylated forms of signaling proteins and kinases, transcription factors, markers of cell state including quiescence, senescence, apoptosis, stress, etc. as well as a variety of non-antibody-based fluorescent stains (Table 2). Multiplexing antibodies and stains makes it possible to discriminate among proliferating, quiescent and dying cells, identify tumor and stroma, and collect immuno-phenotypes (Angelo et al., 2014;Giesen et al., 2014;Goltsev, 2017). Use of phospho-specific antibodies and antibodies against proteins that re-localize upon activation (e.g. transcription factors) makes it possible to assay the states of signal transduction networks. For example, in a 10-cycle t-CyCIF analysis of human tonsil ( Figure 3A) subcellular features such as membrane staining, Ki-67 puncta (Cycle 1), ring-like staining of the nuclear lamina (Cycle 6) and nuclear exclusion of NF-KB (Cycle 6) can easily be demonstrated ( Figure 3B). The five-cycle t-CyCIF data on normal skin in Figure 3C shows tight localization of auto-fluorescence (likely melanin) to the epidermis prior to prebleaching and images of three non-antibody stains used in the last t-CyCIF cycle: HCS CellMask Red Stain for cytoplasm and nuclei, Actin Red, a Phalloidin-based stain for actin and Mito-tracker Green for mitochondria.
anti-S100 S100    In the current work, we rely exclusively on commercial antibodies that have previously been validated using IHC or conventional immunofluorescence; when feasible we confirm that staining by t-CyCIF resembles what has previously been reported for IHC staining. This does not constitute a sufficient level of testing or validation for discovery science or clinical studies and the patterns of staining described in this paper should therefore be considered illustrative of the t-CyCIF approach rather than definitive descriptions; we are currently developing a database of matched t-CyCIF and IHC images across multiple tissues and knockdown cell lines to address this issue and share validation test data with the wider research community.

Fluorophore inactivation, cycle count and tissue integrity
The efficiency of fluorophore inactivation by hydrogen peroxide, light and high pH varies with fluorophore but only minimally with the antibody to which the fluorophore is coupled (Alexa Fluor 488 is inactivated more slowly than Alexa Fluor 570 or 647; Figure 4B and Figure 4-figure supplement 1). We typically incubate specimens in bleaching conditions for 60 min, which is sufficient to reduce fluorescence intensity by 10 2 to 10 3 -fold ( Figure 4C). When testing new antibodies or analyzing new tissues, imaging is performed after each bleaching step and prior to initiation of another t-CyCIF cycle to ensure that fluorophore inactivation is complete. In preliminary studies, we have tested a range of other fluorophores for their compatibility with t-CyCIF including FITC, TRITC, phycoerythrin, Allophycocyanin, eFluor 570 and eFluor 660 (eBioscience). We conclude that it will be feasible to increase the number of t-CyCIF channels per cycle from four to at least six (3 to 5 antibodies plus a DNA stain). However, all the images in this paper are collected using a four-channel method.
The primary limitation on the number of t-CyCIF cycles that can be performed is the integrity of the tissue: some tissues samples are physically more robust and can withstand more staining and washing procedures than others ( Figure 4D). To study the effect of cycle number on tissue integrity, we performed a 10-cycle t-CyCIF experiment on a tissue microarray (TMA) comprising a total of 40 cores from 16 different tissues and tumor types. After each t-CyCIF cycle, the number of nuclei remaining was quantified for each core relative to the initial number. For example, Figure 4D shows breast, bladder, lung and prostate cores in which cell number was reduced after 10 cycles by~2% and an unusually high 46% (apparent increases in cell number in these data are caused by fluctuation in the performance of cell segmentation routines and are not statistically significant). Cells that were lost appear red in these images. The data show that cell loss is often uneven across samples, preferentially affecting regions of tissue with low cellularity.
Overall, we found that the extent of cell loss varied with tissue type and, within a single tissue type, from core to core (six breast cores are shown; Figure 4E). For many tissues, we have not yet   attempted to optimize cycle number and the experiments performed to date do not fully control for pre-analytical variables (Vassilakopoulou et al., 2015) such as fixation time and the age of tissue blocks. As a rule, we find that normal tonsil, skin, glioblastoma, ovarian cancer, pancreatic cancer and melanoma can be subjected to >15 cycles with less than 25% cell loss. Figure 4F shows a melanoma specimen subjected to 20 t-CyCIF cycles with good preservation of cell and tissue morphology ( Figure 4G). We conclude that t-CyCIF is compatible with multiple normal tissues and tumor types but that some tissues and/or specimens can be subjected to more cycles than others. One requirement for high cycle number appears to be cellularity: samples in which cells are very sparse tend to be more fragile. We expect improvements in cycle number with additional experimentation and the use of fluidic devices that deliver staining and wash liquids more gently.
One potential concern about cyclic immunofluorescence is that the process is relatively slow; each cycle takes 6-8 hr and we typically perform one cycle per day. However, a single operator can easily process 30 slides in parallel, and in the case of TMAs, 30 slides can comprise over 2000 different samples. Under these conditions, the most time-consuming step in t-CyCIF is collecting the 200-400 fields of view needed to image each slide. Time could be saved by imaging fewer cells per sample, but the results described below (demonstrating substantial cellular heterogeneity in a single piece of a tumor resection) strongly argue in favor of analyzing as large a fraction of each tissue specimen as possible. As a practical matter, data analysis and data interpretation remain more timeconsuming than data collection. We also note that the throughput of t-CyCIF compares favorably with other tissue-imaging platforms or single-cell transcriptome profiling.

Impact of cycle number on immunogenicity
Because t-CyCIF assembles multiplex images sequentially, it is sensitive to factors that alter immunogenicity as cycle number increases. To investigate such effects, we performed a 16-cycle t-CyCIF experiment in which the order of antibody addition was varied between two immediately adjacent tissue slices cut from the same tissue block ( Figure 5A; Slides A and B); the study was repeated three times, once with tonsil and twice with melanoma specimens with similar results (~1.8 Â 10 5 cells were used for the analysis and overall cell loss was <15%).
This experiment made it possible to judge: (i) the repeatability of staining a single specimen using the same set of antibodies ( Figure 5A, denoted by yellow highlight) (ii) the similarity of staining between slides A and B (blue highlight) and (iii) the effect of swapping the order of antibody addition (cycle number) between slides A and B (blue lines). Comparisons within a single slide were made on a cell-by-cell basis but because slides A and B contain different cells, comparisons between slides were made at the level of intensity distributions (computed on a per-cell basis following segmentation). The repeatability of staining (as measured in cycles 3, 7, 12 and 16) was performed using  Figure 5B). Repeated staining of the same antigen is expected to saturate epitopes, but we reasoned that this effect would be less pronounced the more abundant the antigen. For PCNA, the correlation in staining intensities across four cycles was high (r = 0.95 to 0.99) and somewhat lower in the case of Vimentin and Tubulin (r = 0.80 to 0.95; Figure 6A; a more extensive comparison is shown in Figure 6-figure supplement 1). When we examined the corresponding images, it was readily apparent that Tubulin, and to a lesser extent Vimentin, stained more intensely in later than in earlier t-CyCIF cycles (see intensity distributions in Figure 6A and images in Figure 6B). When images were scaled to equalize the intensity range (by histogram equalization), staining patterns were indistinguishable across all cycles and loss of cells or specific subcellular structures was not obviously a factor ( Figure 6B, left vs right panels and Figure 6C). Thus, for at least a subset of antibodies, staining intensity increases rather than decreases with cycle number whereas background fluorescence falls. As a consequence, dynamic range, defined here as the ratio of the least to the most intense 5% of pixels, frequently increases with cycle number ( Figure 6A and Figure 6-figure supplement 1). These effects were reproducible across slides A and B in all three experiments performed.
When we compared staining between slides A and B for the same antibodies and cycle number, the overlap in intensity distributions was high (>0.85), demonstrating good sample to sample reproducibility (Zhou and Liu, 2012). The overlap remained high for the majority of antibodies even when they were used in different cycles on slides A and B, but for some antibodies, signal intensity clearly increased or decreased with cycle number (Figure 6D; blue and red outlines). In the case of eight antibodies for which the effect of cycle number was greatest (including tubulin, as discussed above), the overlap in intensity distributions was <0.6 as a consequence of both increases and decreases in staining intensity ( Figure 6E). Overall, we found that the repeatability of staining between two biological samples was highest when the antibodies were used in the same cycle on both samples, lower when the antibodies were used in different cycles on the sample, and lowest when both the order and sample were different ( Figure 6F).
The reasons for changes in staining intensity with cycle number are not known, but the fact that the same changes were observed across multiple experiments (for any single antibody) suggests that they arise not from irreproducibility of the t-CyCIF procedure but rather from changes in epitope accessibility. Even in these cases, it appears that it is absolute intensity rather than morphology that is variable. Thus, while changes in staining intensity with cycle number are a concern for a subset of t-CyCIF antibodies, it should be possible to minimize the problem by staining all samples in the Figure 6 continued values were integrated across whole cells to construct the distribution. Box plots to right: estimated dynamic range at four cycle numbers 3, 7, 12, 16. Red lines denote median intensity values (across 56 frames), boxes denote the upper and lower quartiles, whiskers indicate values outside the upper/ lower quartile within 1.5 standard deviations, and red dots represent outliers. (B) Representative images showing anti-tubulin Alexa 647 staining at four t-CyCIF cycles; original images are shown on the left (representing the same exposure time and approximately the same illumination) and images scaled by histogram equalization to similar intensity ranges are shown on the right. (C) Image for anti-CD45RO-Alexa 555 at cycles 5 and 15 scaled to similar intensity ranges as described in (B); the dynamic range (DR) of the cycle 15 image is~3.3 fold lower than that of the Cycle 5 image, but shows similar morphology. (D) Intensity distributions for selected antibodies that were used in different cycles on Slides A and B. Colors denote the degree of concordance between the slides ranging from high (overlap >0.8 in yellow; PCNA), slightly increased or decreased with increasing cycle (overlap 0.6 to 0.8 in light blue or light red; S100 and SMA) or substantially increased or decreased (overlap <0.6 in red or blue; VEGFR2 and CD45RO). (E) Summary of effects of cycle number on antibody staining based on the degree of overlap in intensity distributions (the overlap integral); color coding is the same as in (D). (F) Effect of cycle number and specimen identity on overlap integrals for all antibodies and all cycles assayed. The red line denotes the median intensity value, boxes denote the upper/lower quartiles, and whiskers indicate values outside the upper/lower quartile and within 1.5 standard deviations, and red dots represent outliers. All the numeric data in Figures 5 and 6 Figure 7 continued on next page same order. Other approaches will also be important; for example, using calibration standards and identifying antibodies exhibiting the least variation with cycle number.
One way to reduce artefacts generated by differences in the order of antibody addition is to create a single high-plex antibody mixture and then stain all antigens in parallel. This approach is not compatible with t-CyCIF but is feasible using methods such as MIBI or CODEX (Angelo et al., 2014;Goltsev, 2017). However, there is substantial literature showing that the formulation of highly multiplex immuno-assays is complicated by interaction among antibodies (Ellington et al., 2010) that has a physicochemical explanation in some cases in weak self-association and viscosity . Consistent with these data, we have observed that when eight or more unlabeled antibodies are added to a t-CyCIF experiment, the intensity of staining can fall, although the effect is smaller than observed with antibodies most sensitive to order of addition. We conclude that the construction of sequentially applied t-CyCIF antibody panels and of single high-plex mixtures will both require optimization of specific panels and their method of use.

Analysis of large specimens by t-CyCIF
Review of large histopathology specimens by pathologists involves rapid and seamless switching between low-power fields to scan across large regions of tissue and high-power fields to study cellular morphology. To mimic this integration of information at both tissue and cellular scales, we performed eight-cycle t-CyCIF on a large 2 Â 1.5 cm resection specimen that includes pancreatic ductal adenocarcinoma (PDAC) and adjacent normal pancreatic tissue and small intestine ( Figure 7A-C). Nuclei were located in the DAPI channel and cell segmentation performed using a watershed algorithm (Figure 7-figure supplement 1: see Materials and methods section for a discussion of the method and its caveats) yielding~2 Â 10 5 single cells each associated with a vector comprising 25 whole-cell fluorescence intensities. Differences in subcellular distribution were evident for many proteins, but for simplicity, we only analyzed fluorescence intensity on a per-antigen basis integrated over each whole cell. Results were visualized by plotting intensity value onto the segmentation data ( Figure 7D), by computing correlations on a cell-by-cell basis ( Figure 7E), or by using t-distributed stochastic neighbor embedding (t-SNE) (Maaten and Hinton, 2008), which clusters cells in 2D based on their proximity in the 25-dimensional space of image intensity data ( Figure 8A).
The analysis in Figure 7E shows that E-cadherin, keratin and b-catenin levels are highly correlated with each other, whereas vimentin and VEGFR2 receptor levels are anti-correlated, recapitulating the known dichotomy between epithelial and mesenchymal cell states in normal and diseased tissues. Many other physiologically relevant correlations are also observed, for example between the levels of pERK T202/Y204 (the phosphorylated, active form of the kinase) and activating phosphorylation of the downstream kinase pS6 S235/S236 (r = 0.81). When t-SNE was applied to all cells in the specimen, we found that those identified during histopathology review as being from non-neoplastic pancreas (red) were distinct from PDAC (green) and also from the neighboring non-neoplastic small intestine (blue) ( Figure 8B-D). Vimentin and E-Cadherin had very different levels of expression in PDAC and normal pancreas as a consequence of epithelial-to-mesenchymal transitions (EMT) in malignant tissues as well as the presence of a dense tumor stroma, a desmoplastic reaction that is a hallmark of the PDAC microenvironment (Mahadevan and Von Hoff, 2007). The microenvironment  Figure 8. DOI: https://doi.org/10.7554/eLife.31657.021

(E) Quantitative single-cell signal intensities of 24 proteins (rows) measured in~4Â10 3 cells (columns) from panel (C). The Pearson correlation coefficient for each measured protein with E-cadherin (at a single-cell level) is shown numerically. Known dichotomies are evident such as anti-correlated expression of epithelial (E-Cadherin) and mesenchymal (Vimentin) proteins. Proteins highlighted in red are further analyzed in
The following source data and figure supplement are available for figure 7: Source data 1. Single-cell intensity data used in Figure 7E.  Figure 7 with the fluorescence intensities for markers of proliferation (PCNA and Ki67) and signaling (pERK and b-catenin) overlaid on the plots as heat maps. In both tissue types, there exists substantial heterogeneity: circled areas indicate the relationship between pERK and b-catenin levels in cells and represent positive ('a'), negative ('b') or no association ('c') between these markers. (B) Figure 8 continued on next page of PDAC was more heavily infiltrated with CD45 + immune cells than the normal pancreas, and the intestinal mucosa of the small intestine was also replete with immune cells, consistent with the known architecture and organization of this tissue.
The capacity to image samples that are several square centimeters in area with t-CyCIF can facilitate the detection of signaling biomarker heterogeneity. The WNT pathway is frequently activated in PDAC and is important for oncogenic transformation of gastrointestinal tumours (Jones et al., 2008). Approximately 90% of sporadic PDACs also harbor driver mutations in KRAS, activating the MAPK pathway and promoting tumourigenesis (Vogelstein et al., 2013). Studies comparing these pathways have come to different conclusions with respect to their relationship: some studies show concordant activation of MAPK and WNT signaling and others argue for exclusive activation of one pathway or the other (Jeong et al., 2012). In t-SNE plots derived from images of PDAC, multiple sub-populations of cells representing negative, positive or no correlation between pERK and b-catenin levels can be seen (marked with labels 'a', 'b' or 'c', respectively in Figure 8A). The same three relationships can be found in non-neoplastic pancreas and small intestine ( Figures 8A and 7C). In PDAC, malignant cells can be distinguished from stromal cells, to a first approximation, by high proliferative index, which can be measured by staining for Ki-67 and PCNA (Bologna-Molina et al., 2013). When we gated for cells that were both Ki67 high and PCNA high , and thus likely to be malignant, the co-occurrence of different relationship between pERK and b-catenin levels on a cellular level was again evident. While we cannot exclude the possibility of phospho-epitope loss during sample preparation, it appears that the full range of possible relationships between the MAPK and WNT signaling pathways described in the literature can be found within a specimen from a single patient, illustrating the impact of tissue context on the activities of key signal transduction pathways.

Multiplex imaging of immune infiltration
Immuno-oncology drugs, including immune checkpoint inhibitors targeting CTLA-4 and the PD-1/ PD-L1 axis are rapidly changing the therapeutic possibilities for traditionally difficult-to-treat cancers including melanoma, renal and lung cancers, but responses are variable across and within cancer types. The hope is that tumor immuno-profiling will yield biomarkers predictive of therapeutic response in individual patients. For example, expression of PD-L1 correlates with responsiveness to the ICIs pembrolizumab and nivolumab (Mahoney and Atkins, 2014) but the negative predictive value of PD-L1 expression alone is insufficient to stratify patient populations (Sharma and Allison, 2015). In contrast, by measuring PD-1, PD-L1, CD4 and CD8 by IHC on sequential tumor slices, it has been possible to identify some immune checkpoint inhibitor-responsive melanom patients (Tumeh et al., 2014). To test t-CyCIF in this application, eight-cycle imaging was performed on a 1 Â 2 cm specimen of clear-cell renal cell carcinoma using 10 antibodies against multiple immune markers and 12 against other proteins expressed in tumor and stromal cells ( Figure 9A-B; Supplementary file 4). A region of the specimen corresponding to tumor was readily distinguishable from non-malignant stroma based on a-SMA expression (a-SMA high regions denote stroma and a-SMA low regions high density of malignant cells).
In the a-SMA low domain, CD3 + or CD8 + lymphocytes were fourfold enriched ( Figure 9C) and PD-1 and PD-L1-positive cells were 13 to 20-fold more prevalent as compared to the surrounding tumor stroma (a-SMA high domain); CD3 + CD8 + double positive T-cells were found almost exclusively in the tumor. Suppression of immune cells is mediated by binding of PD-L1 ligand, which is commonly expressed by tumor cells, to the PD1 receptor expressed on immune cells (Tumeh et al., 2014). To Quantitative assessment of total lymphocytic cell infiltrates (CD3 + cells), CD8 + T lymphocytes, cells expressing PD-1 or its ligand PD-L1 or the VEGFR2 Figure 9 continued on next page begin to estimate the likelihood of ligand-receptor interactions, we quantified the degree of colocalization of cells expressing the two molecules. The centroids of PD-1 + or PD-L1 + cells were determined from images (PD-1, red; PD-L1, green, Figure 9E) and co-localization (highlighted in yellow, Figure 9F) computed by k-nearest neighbor analysis. We found that co-localization of PD-1/PD-L1 was~2.7-fold more likely (Figure 9-figure supplement 1) in tumor and stroma and was concentrated on the tumor-stroma border consistent with previous reports on melanoma (Tumeh et al., 2014). These data demonstrate the potential of spatially resolved immuno-phenotyping to quantify state and location of tumor infiltrating lymphocytes; such data may ultimately yield biomarkers predictive of sensitivity to immune checkpoint inhibitor (Tumeh et al., 2014).

Analysis of diverse tumor types and grades using t-CyCIF of tissuemicroarrays (TMA)
To explore the general utility of t-CyCIF in a range of healthy and cancer tissues we applied eight cycle t-CyCIF to TMAs containing 39 different biopsies from 13 healthy tissues and 26 biopsies corresponding to low-and high-grade cancers from the same tissue types ( Figure 10A and Figure 10figure supplement 1, Supplementary file 3 for antibodies used, Supplementary file 5 for TMA details and naming conventions) and then performed t-SNE and clustering on single-cell intensity data ( Figure 10B). The great majority of TMA samples mapped to one or a few discrete locations in the t-SNE projection (compare normal kidney tissue -KI1, low-grade tumors -KI2, and high-grade tumors -KI3; Figure 10C), although ovarian cancers were scattered across the t-SNE projection ( Figure 10D); overall, there was no separation between normal tissue and tumors regardless of grade ( Figure 10E). In a number of cases, high-grade cancers from multiple different tissues of origin co-clustered, implying that transformed morphologies and cell states were closely related. For example, while healthy and low-grade pancreatic and stomach cancer occupied distinct t-SNE domains, high-grade pancreatic and stomach cancers were intermingled and could not be readily distinguished ( Figure 10F), recapitulating the known difficulty in distinguishing high-grade gastrointestinal tumors of diverse origin by histophathology (Varadhachary and Raber, 2014). Nonetheless, t-CyCIF might represent a means to identify discriminating biomarkers by efficiently sorting through large numbers of alternative antigens and antigen localizations.
Quantitative analysis reveals global and regional heterogeneity and multiple histologic subtypes within the same tumor in glioblastoma multiforme (GBM) Data from single-cell genomics reveals extensive heterogeneity in many types of cancer (Turner and Reis-Filho, 2012) but our understanding of this phenomenon requires spatially resolved data (Giesen et al., 2014). We performed eight-cycle imaging on a 2.5 cm x 1.8 mm resected glioblastoma (GBM) specimen imaging markers of neural development, cell cycle state and signal transduction ( Figure 11A-B, Supplementary file 6). GBM is a highly aggressive and genetically heterogeneous (Brennan et al., 2013) brain cancer commonly classified into four histologic subtypes Figure 9 continued for the entire tumor or for a-SMA high and a-SMA low regions. VEGFR2 is a protein primarily expressed in endothelial cells and is targeted in the treatment of renal cell cancer. The error bars represent the S.E.M. derived from 100 rounds of bootstrapping. (D) Density plot for CD3 and CD8 expression on single cells in the tumor (left) or stromal domains (right). (E) Centroids of CD3 + or CD3 + CD8 + cells in blue or dark blue as well as cells staining as SMA high or SMA low (gray and light-gray, respectively) used to define the stromal and tumor regions. (F) Centroids of PD-1 + and PD-L1 + cells are shown in red and green, respectively. (G) Results of a K-nearest neighbor algorithm used to compute areas in which PD-1 + and PD-L1 + cells lie within~10 mm of each other and with high spatial density (in yellow) and thus, are potentially positioned to interact at a molecular level. DOI: https://doi.org/10.7554/eLife.31657.027 The following source data and figure supplement are available for figure 9: Source data 1. Immune cell counts from bootstrapping in tumor and stroma regions ( Figure 9C). DOI: https://doi.org/10.7554/eLife.31657.029 Source data 2. Single-cell intensity data used in Figure 9.  (Olar and Aldape, 2014). Following image segmentation, phenotypic heterogeneity was assessed at three spatial scales corresponding to: (i) 1.6 Â 1.4 mm fields of view (252 total) each of which comprised 10 3 to 10 4 cells (ii) seven macroscopic regions of~10 4 to 10 5 cells each, corresponding roughly to tumor lobes and (iii) the whole tumor comprising~10 6 cells. To quantify local heterogeneity, we computed the informational entropy on a-per-channel basis for 10 3 randomly selected cells in each field ( Figure 11C; see online Materials and methods for details). In this setting, informational entropy is a measure of cell-to-cell heterogeneity on a mesoscale corresponding to 10-30 cell diameters. For a marker such as EGFR, which can function as a driving oncogene in GBM, informational entropy was high in some areas ( Figure 11C; red dots) and low in others (blue dots). Areas with high entropy in EGFR abundance did not co-correlate with areas that were most variable with respect to a downstream signaling protein such as pERK. Thus, the extent of local heterogeneity varied with the region of the tumor and the marker being assayed. Semi-supervised clustering using expectation-maximization Gaussian mixture (EMGM) modeling of all cells in the tumor yielded eight distinct clusters, four of which encompassed 85% of all cells ( Figure 12A and Figure 12-figure supplement 1). Among these, cluster one had high EGFR levels, cluster two had high NGFR and Ki67 levels and cluster six had high levels of vimentin; cluster five was characterized by high keratin and pERK levels. The presence of four highly populated t-CyCIF clusters is consistent with data from single-cell RNA-sequencing of~400 cells from five GBMs (Patel et al., 2014). Three of the t-CyCIF clusters have properties reminiscent of established histological subtypes including: classical, cluster 1; pro-neural, cluster 3; and mesenchymal, cluster 6, but additional work will be required to confirm such assignments.
To study the relationship between phenotypic diversity and tumor architecture, we mapped each cell to an EMGM cluster (denoted by color). Extensive intermixing was observed at all spatial scales ( Figure 12B). For example, field of view 147 was highly enriched for cells corresponding to cluster 5 (yellow), but a higher magnification view revealed extensive intermixing of four other cluster types on a scale of~3-5 cell diameters ( Figure 12C). At the level of larger, macroscopic tumor regions, the fraction of cells from each cluster also varied dramatically ( Figure 12D). None of these findings was substantially different when the number of clusters was set to 12 ( Figure 12-figure supplement 2).
These results have several implications. First, they suggest that GBM is phenotypically heterogeneous on a spatial scale of 5-1000 cell diameters and that cells corresponding to distinct t-CyCIF clusters are often found in the vicinity of each other. Second, sampling a small region of a large tumor has the potential to misrepresent the proportion and distribution of tumor subtypes, with implications for prognosis and therapy. Similar concepts likely apply to other tumor types with high genetic heterogeneity, such as metastatic melanoma (Tirosh et al., 2016), and are therefore relevant to diagnostic and therapeutic challenges arising from tumor heterogeneity. Figure 10 continued data were analyzed using the CYT package (see Materials and methods). Tissues of origin and corresponding malignant lesions were labeled as follows: BL, bladder cancer; BR, breast cancer CO, Colorectal adenocarcinoma, KI, clear cell renal cancer, LI, hepatocellular carcinoma, LU, lung adenocarcinoma, LY, lymphoma, OV, high-grade serous adenocarcinoma of the ovary, PA, pancreatic ductal adenocarcinoma, PR, prostate adenocarcinoma, UT, uterine cancer, SK, skin cancer (melanoma), ST, stomach (gastric) cancer. Numbers refer to sample type; '1' to normal tissue, '2' to -grade tumors and '3' to high-grade tumors. (C) Detail from panel B of normal kidney tissue (KI1) a low-grade tumor (KI2) and a high-grade tumor (KI3) (D) Detail from panel B of normal ovary (OV1) low-grade tumor (OV2) and high-grade tumor (OV3). (E) t-SNE plot from Panel B coded to show the distributions of all normal, low-grade and high-grade tumors. (F) tSNE clustering of normal pancreas (PA1) and pancreatic cancers (low-grade, PA2, and high-grade, PA3) and normal stomach (ST1) and gastric cancers (ST2 and ST3, respectively) showing intermingling of high-grade cells. DOI: https://doi.org/10.7554/eLife.31657.031 The following source data and figure supplement are available for figure 10: Source data 1. Single-cell intensity data used in Figure 10. DOI: https://doi.org/10.7554/eLife.31657.033 Figure supplement 1. Gallery of exemplary tissues imaged on the TMA described in Figure 10. DOI: https://doi.org/10.7554/eLife.31657.032

Discussion
The complex molecular biology and spatial organization of tissues and solid tumors poses a scientific and diagnostic challenge that is not sufficiently addressed using single-cell genomics, in which morphology is commonly lost, or H&E and single-channel IHC staining, which provide data on only a few Figure 11. Molecular heterogeneity in a single GBM tumor. (A) Representative low-magnification image of a GBM specimen generated from 221 stitched 10X frames; the sample was subjected to 10 rounds of t-CyCIF using antibodies listed in Supplementary file 6. (B) Magnification of frame 152 (whose position is marked with a white box in panel A) showing staining of pERK, pRB and EGFR; lower panel shows a further magnification to allow single cells to be identified. (C) Normalized Shannon entropy of each of 221 fields of view to determine the extent of variability in signal intensity for 1000 cells randomly selected from that field for each of the antibodies shown. The size of the circles denotes the number of cells in the field and the color represents the value of the normalized Shannon entropy (data are shown only for those fields with more than 1000 cells; see Materials and methods for details). DOI: https://doi.org/10.7554/eLife.31657.034 The following source data is available for figure 11: Source data 1. Normalized entropy data shown in Figure 11C. DOI: https://doi.org/10.7554/eLife.31657.035 Source data 2. Single-cell intensity data used in Figure 11 and 12. DOI: https://doi.org/10.7554/eLife.31657.036 proteins or molecular features. At the same time, the vast number of FFPE histological specimens collected in the course of routine clinical care and clinical trials (and in the study of model organisms) represents an underutilized resource with great potential for novel discovery. A variety of methods for performing highly multiplexed immune-based imaging of cells and tissues has recently been described including imaging cytometry (Giesen et al., 2014), MIBI (Angelo et al., 2014), DNAexchange imaging (DEI) (Wang, 2017) and CODEX (Goltsev, 2017); FISSEQ (Lee et al., 2014) directly images expressed RNAs. Like traditional antibody stripping approaches, the cyclic immunofluorescence approach first described by Gerdes et al (Gerdes et al., 2013) and further developed here assembles highly multiplexed images by sequential acquisition of lower dimensional immunofluorescence images. We show here that the t-CyCIF implementation of cyclic immunofluorescence is compatible with a wide range of antibodies and tissue types and yields up to 60-plex images with excellent preservation of small intracellular structures.
The requirement in t-CyCIF for multiple rounds of staining and imaging might seem to be a liability but it has several substantial advantages relative to all-in-one methods such as MIBI, DEI and CODEX. First, t-CyCIF can be performed using existing fluorescence microscopes. Not only does this reduce costs and barriers to entry, it allows the unique strengths of slide-scanning, confocal, and structured illumination microscopes to be exploited. Using different instruments, samples several square centimeters in area can be rapidly analyzed at resolutions of~1 mm and selected fields of view studied at super-resolution (~110 nm on an OMX Blaze). Multiscale imaging makes it possible to combine tissue-level architecture with subcellular morphology, much like a pathologist switching between low-and high-power fields, but there is little chance that such capabilities can be combined in a single instrument. Because no spectral deconvolution is required, t-CyCIF can use highly optimized filter sets and fluorophores, resulting in good sensitivity. t-CyCIF antibody panels are also simple to assemble and validate using commercial antibodies, including those that constitute FDAapproved diagnostics. This avoids the limitations of an exlusive reliance on pre-assembled reagent kits provided by manufacturers. Finally, t-CyCIF is compatible with H&E staining, enabling fluorescence imaging to be combined with conventional histopathology review.
Commercial systems for non-optical tissue imaging are only now starting to appear and it is difficult to compare their performance to multiplexed immunofluorescence, particularly because the approach published by Gerdes et al. (2013) is proprietary and available only as commercial service. In contrast, the t-CyCIF method described here can easily be implemented in a conventional research or clinical laboratory without the need for expensive equipment or specialized reagents. As MIBI, DEI and CODEX instruments come on-line, direct comparison with t-CyCIF will be possible. We anticipate that high resolution and good linearity will be areas in which fluorescence imaging is superior to enzymatic amplification, laser ablation or mechanical picking of tissues. t-CyCIF is relatively slow when performed on a single sample, but when many large specimens or TMAs are processed in parallel, throughput is limited primarily by imaging acquisition, which is at least as fast as approaches involving laser ablation. Considerable opportunity exists for further improvement in t-CyCIF by switching from four to six-channels per cycle, optimizing bleach and processing solutions to preserve tissue integrity, using fluidic devices to rapidly process many slides in parallel and developing better software for identifying fields of view that can be skipped in large irregular specimens. Because direct fluorescence will remain challenging in the case of very rare epitopes, we speculate that hybrid approaches involving t-CyCIF and methods such as DEI or CODEX will ultimately prove to be most effective.
As in all methods involving immune detection, antibodies are the most critical and difficult to validate reagents in t-CyCIF. To date, we have shown that over 200 commercial antibodies are compatible with the method as judged by patterns of staining similar to those previously reported for IHC; this is an insufficient level of validation for most studies and we are therefore working to develop a generally useful antibody validation resource (www.cycif.org). Thus, while this paper describes markers relevant to diagnosis of disease, our results are illustrative of the t-CyCIF approach and specific findings might not prove statistically significant when tested on larger, well-controlled sets of human samples.
There is little or no evidence that antigenicity falls across the board in t-CyCIF as cycle number increases; signal-to-noise ratios can even increase due to falling background auto-fluorescence. When samples are stained with the same antibodies in different t-CyCIF cycles, repeatability is high (as measured by correlation in staining intensity on a cell-by-cell basis) as is reproducibility across  Figure 11. Intensity values were clustered using expected-maximization with Gaussian mixtures (EMGM), yielding eight clusters, of which four clusters accounted for the majority of cells. The intensity scale shows the average level for each intensity feature in that cluster. The number of cells in the cluster is shown as a percentage of all cells in the tumor (bottom of panel). An analogous analysis is shown for 12 clusters in Figure 12 two successive slices of tissue (as measured by overlap in intensity distributions). Moreover, for the majority of antibodies tested, order of use is not critical. For some antibodies fluorescence intensity increases with cycle number and for others it decreases; these factors need to be considered when developing a staining strategy. While the precise reasons for variation in staining with cycle number are not known such variation is reproducible across specimens, suggesting that it reflects properties of the epitope or antibody and not the t-CyCIF process per se , variation in staining can be minimized by staining all specimens with the same antibodies in the same order (which also represents the most practical approach). However, this solution is likely to be insufficient for creation of largescale t-CyCIF datasets in which diverse tissues will be compared with each other (e.g. in proposed tissue atlases [Department of Health and Human Services, 2018]) and it will therefore be important to identify antibodies for which cycle number has minimal impact and to create effective methods to correct for those fluctuations that do occur (e.g. inclusion of staining controls).
As an initial application of t-CyCIF, we examined a cancer resection specimen that includes PDAC, healthy pancreas and small intestine. Images were segmented and fluorescence intensities in~10 5 whole cells calculated for 24 antibody channels plus a DNA stain. Integrating intensities in this manner does not make use of the many subcellular features visible in t-CyCIF images and therefore represents only a first step in data analysis. We find that expression of vimentin and E-cadherin, classical markers of epithelial and mesenchymal cells, are strongly anti-correlated at a single-cell level and that malignant tissue is skewed toward EMT, consistent with prior knowledge on the biology of pancreatic cancer (Zeitouni et al., 2016). The WNT and ERK/MAPK pathways are known to play important roles in the development of PDAC (Jones et al., 2008), but the relationship between the two pathways remains controversial. t-CyCIF reveals a negative correlation between b-catenin levels (a measured of WNT pathway activity) and pERK (a measure of MAPK activity) in cells found in some regions of PDAC, non-malignant small intestine and pancreas, a positive correlation in other regions and no significant correlation in yet others. Thus, the full range of discordant observations found in the literature can be recapitulated within a single tumor, emphasizing the wide diversity of signaling states observable at a single-cell level.
As a second application of t-CyCIF, we studied within-tumor heterogeneity in GBM, a brain cancer with multiple histological subtypes whose differing properties impact prognosis and therapy (Olar and Aldape, 2014;Phillips et al., 2006). Clustering reveals multiple phenotypic classes intermingled at multiple spatial scales with no evidence of recurrent patterns. In the GBM we have studied in detail, heterogeneity on a scale of 10-100 cell diameters is as great as it is between distinct lobes. The proportion of cells from different clusters also varies dramatically from one tumor lobe to the next. Although it is not yet possible to link t-CyCIF clusters and known histological subtypes, cell-to-cell heterogeneity on these spatial scales are likely to impact the interpretation of small biopsies (e.g. a core needle biopsy) of a large tumor sample; the data also emphasize the inherent limitation in examining only a small part of a large tumor specimen (e.g. to save time on image acquisition). At the same time, it is important to note that cell-to-cell heterogeneity is caused by processes operating on a variety of time scales, only some of which are likely to be relevant to therapeutic response and disease progression. For example, some cell-to-cell differences visible in GBM images arise from a cyclic process, such as cell cycle progression, whereas others appear to involve differences in cell lineage or clonality. Methods to correct for the effects of variation in cell cycle state have been worked out for single-cell RNA-sequencing (Izar, 2017), but will require further work in imaging space. In a third application of t-CyCIF, we characterized tumor-immune cell interactions in a renal cell tumor. Immune checkpoint inhibitors elicit durable responses in a portion of patients with diverse types of cancer, but identifying potential responders and non-responders remains a challenge. In those cancers in which it has been studied (Mahoney and Atkins, 2014), quantification of single checkpoint receptors or ligands by IHC lacks sufficient positive and negative predictive value to stratify therapy or justify withholding checkpoint inhibitors in favor of small molecule therapy (Sharma and Allison, 2015). Multivariate predictors based on multiple markers such as CD3, CD4, CD8, PD-1 etc. appear to be more effective, but still underperform in patient stratification (Tumeh et al., 2014) probably because cells other than CD8 +lymphocytes affect therapeutic responsiveness. In this paper, we perform a simple analysis to show that tumor infiltrating lymphocytes can be subtyped by t-CyCIF and analyzed for the proximity of PD-1 and PD-L1 at a single-cell level. Next steps involve thorough interrogation of immuno-phenotypes by multiplex imaging to relate staining patterns in images to immune cell classes previously defined by flow cytometry and to identify immune cell states that fall below the limit of detection for existing analytical methods.
In conclusion, t-CyCIF is a robust, easy to implement approach to multi-parametric tissue imaging applicable to many types of tumors and tissues; it allows investigators to mix and match antibodies Table 3. Breakdown of individual steps performed for dewaxing and antigen retrieval on a Leica BOND.
Step Reagent Supplier Incubation (min) Temp. (˚C) depending on the requirements of a specific type of sample. To create a widely available community resource, we have posted antibody lists, protocols and example data at http//www.cycif.org and are currently updating this information on a regular basis. Highly multiplexed histology is still in an early stage of development and better methods for segmenting cells, quantifying fluorescence intensities and analyzing the resulting data are in development by multiple groups. The resulting ability to quantify cell-to-cell heterogeneity may enable reconstruction of signaling network topologies in situ (Giesen et al., 2014;Sachs et al., 2002) by exploiting the fact that protein abundance and states of activity fluctuate from one cell to the next; when fluctuations are well correlated, they are likely to reflect causal associations (Vilela and Danuser, 2011). We expect t-CyCIF to be complementary to, and used in parallel with other protein and RNA imaging methods such as FISSEQ (Lee et al., 2015) or DEI  that may have higher sensitivity or greater channel capacity. A particularly important task will be cross-referencing tumor cell types identified by singlecell genomics or multi-color flow cytometry with those identified by multiplexed imaging, making it possible to precisely define the genetic geography of human cancer and infiltrating immune cells.

Manual dewaxing, rehydration and pre-staining
In our experience dewaxing, rehydration and pre-staining can also be performed manually with similar results. For manual pre-processing, FFPE slides were first incubated in a 60˚C oven for 30 min. To completely remove paraffin, slides were placed in a glass slide rack and then immediately immersed in Xylene in a glass staining dish (Wheaton 900200) for 5 min and subsequently transferred to another dish containing fresh Xylene for 5 min. Rehydration was achieved by sequentially immersing slides, for 3 min each, in staining dishes containing 100% ethanol, 90% ethanol, 70% ethanol, 50% ethanol, 30% ethanol, and then in two successive 1xPBS solutions. Following rehydration, slides were placed in a 1000 ml beaker filled with 500 ml citric acid, pH 6.0, for antigen retrieval. The beaker containing slides and citric acid buffer was microwaved at low power until the solution was at a boiling point and maintained at that temperature for 10 min. After cooling to room temperature, slides were washed 3 times with 1xPBS in vertical staining jars.

Prestaining
Dewaxed specimens were blocked by incubation with Odyssey blocking buffer for 30 mins by applying the buffer to slides as a 250-500 ml droplet at room temperature; evaporation was minimized by using a slide moisture chamber (Scientific Device Laboratory,. Slides were then pre-stained by incubation with diluted secondary antibodies (listed above) for 60 min, followed by washing three times with 1xPBS. Finally, slides were incubated with Hoechst 33342 (2 mg/ml) in 250-500 ml Odyssey blocking buffer for 30 min in a moisture chamber and washed three times with 1xPBS in vertical staining jars. After imaging, cells were subjected to a round of fluorophore inactivation (see below). Following fluorophore inactivation, slides were washed four times with 1x PBS by dipping them in a series of vertical staining jars to remove residual inactivation solution.

Performing cyclic immunofluorescence
All primary antibodies (fluorophore-conjugated and unconjugated) were diluted in Odyssey blocking buffer. Slides carrying tissues that had been subjected to pre-staining, or to a previous t-CyCIF stain and bleach cycle, were incubated at 4˚C for~12 hr with diluted primary or fluorophore-conjugated antibody (250-500 ml per slide) in a moisture chamber. Long incubation times were a matter of convenience and many antibodies only require short incubation with sample. Slides were then washed four times in 1x PBS by dipping in a series of vertical staining jars. For indirect immunofluorescence, slides were incubated in diluted secondary antibodies in a moisture chamber for 1 hr at room temperature followed by four washes with 1xPBS. Slides were incubated in Hoechst 33342 at 2 mg/ml in Odyssey blocking buffer for 15 min at room temperature, followed by four washes in 1xPBS. Stained slides were mounted prior to image acquisition (see the Mounting section below).

Primary antibodies
For t-CyCIF, we selected commercial antibodies previously validated by their manufacturers for use in immunofluorescence, immunocytochemistry or immunohistochemistry (IF, ICC or IHC). When possible, we checked antibodies on reference tissue known to express the target antigen, such as immune cells in tonsil tissue or tumor-specific markers in tissue microarrays. The staining patterns for antibodies with favorable signal-to-noise ratios were compared to those previously reported for that antigen by conventional antibodies. An updated list of all antibodies tested to date can be found at http://www.cycif.org. In current practice, the degree of validation is quantified on a level between 0 and 2: 'Level 0' represents antibodies with inconsistent or no staining in tissues for which the antigen is thought to be present based on published data; 'Level 1' represents the expected pattern of positive staining in a limited number of tissues types (e.g. CD4 antibody in tonsil tissue alone); 'Level 2' represents the expected pattern of positive staining in all tissues or tumor types tested (N >= 3). Higher levels will be assigned in the future to antibodies that have undergone extensive validation; for example, side-by-side comparison of against an established IHC positive control. Overall, the validation of primary antibodies used in this study is not meaningfully greater what has already been done by commercial vendors using conventional IF or IHC.

Mounting and de-coverslipping
Immediately prior to imaging, slides were mounted with 1xPBS or, if imaging was expected to take longer than 30 min, for example, in the case of samples larger than 2-4 cm 2 (corresponding to about 200 fields of view with a 10X objective) PBS was supplement with 10% Glycerol. Slides were covered using 24 Â 60 mm No. one coverslips (VWR 48393-106) to prevent evaporation while facilitating subsequent de-coverslipping via gravity. Following image acquisition, slides were placed in a vertical staining jar containing 1xPBS for at least 15 min. Coverslips were released from slides (and the tissue sample) via gravity as the slides were slowly drawn out of the staining jar.

Fluorophore inactivation (bleaching)
After imaging, fluorophores were inactivated by placing slides horizontally in 4.5% H 2 O 2 and 24 mM NaOH made up in PBS for 1 hr at RT in the presence of white light. Following fluorophore inactivation, slides were washed four times with 1x PBS by dipping them in a series of vertical staining jars to remove residual inactivation solution.

Image acquisition
Stained slides from each round of CyCIF were imaged with a CyteFinder slide scanning fluorescence microscope (RareCyte Inc. Seattle WA) using either a 10X (NA = 0.3) or 40X long-working distance objective (NA = 0.6). Imager5 software (RareCyte Inc.) was used to sequentially scan the region of interest in four fluorescence channels. These channels are referred to by the manufacturer as a: (i) 'DAPI channel' with an excitation filter having a peak of 390 nm and half-width of 18 nm and an emission filter with a peak of 435 nm and half-width of 48 nm; (ii) 'FITC channel' having a 475/28 nm excitation filter and 525/48 nm emission filter (iii); 'Cy3 channel' having a 542/27 nm excitation filter and 597/45 nm emission filter and (iv); 'Cy5 channel' having a 632/22 nm excitation filter and 679/34 nm emission filter. Imaging was performed with 2 Â 2 binning to increase sensitivity, shorten exposure time and reduce photo bleaching. We have tested slide scanners from several other manufacturers (e.g. a Leica Aperio Digital Pathology Slide Scanner, GE IN-Cell Analyzer 6000 and GE Cytell Cell Imaging System) and found that they too can be used to acquire images from samples processed by t-CyCIF. Slides can also be analyzed on conventional microscopes, but the field of view is typically smaller, and an automated stage is required for accurate stitching of individual fields of view into a complete image of a tissue.

Super-resolution microscopy
We acquired 3D-SIM images on a Deltavision OMX V4 Blaze (GE Healthcare) with a 60x/1.42N.A. Plan Apo oil immersion objective lens (Olympus) and three Edge 5.5 sCMOS cameras (PCO). Two to three micron z-stacks were collected with a z-step of 125 nm or 250 nm and with 15 raw images per plane. To minimize spherical aberration, immersion oil matching was used for each sample as described by Hiraoka et al. (1990). except that we measured point spread functions of point-like structures within the sample as opposed to beads on a separate slide. DAPI fluorescence was excited with a 405 nm laser and collected with a 477/35 emission filter, Alexafluor 488 with a 488 nm laser and a 528/48 emission filter, Alexa fluor 555 with a 568 nm laser and a 609/37 emission filter, and Alexa fluor 647with a 642 nm laser and a 683/40 emission filter. All stage positions were saved in softWorX to be revisited later. Super-resolution images were computationally reconstructed from the raw data sets with a channel-specific, measured optical transfer function and a Wiener filter constant of 0.001 using CUDA-accelerated 3D-SIM reconstruction code based on Gustafsson et al. (2008). A comparison of properties of different imaging platforms used in this study are shown in Table 1.

Image processing
Quantitative analysis of tissue images is challenging, in large part because cells are close together and embedded in a complex extracellular environment. Background can be uneven across large images and signal-to-noise ratios relatively low, particularly in the case of tissues with high auto-fluorescence and low signal antibodies (e.g. phospho-protein antibodies). We have only started to tackle these issues in the case of high-dimensional t-CyCIF data and users are encouraged to check for updates on www.cycif.org and implement their own approaches.

Background subtraction and image registration
Background subtraction was performed using the previously established rolling ball algorithm (with a 50-pixel radius) in ImageJ. Adjacent background-subtracted images from the same sample were then registered to each using an ImageJ script as described previously (Lin et al., 2015). All images with 2Â2 binning in acquisition were partially de-convoluted with unsharp masking. DAPI images from each cycle were used to generate reference coordinates by Rigid-body transformation. To generate virtual hyper-stacked images, the transformed coordinates were applied to images from four channel imaging of each t-CyCIF cycle.

Single-cell segmentation and quantification
To obtain intensity values for single cells, images were segmented using a previously described (Lin et al., 2015) Watershed algorithm based on nuclear staining by Hoechst 33342. Images were initially thresholded using the OTSU algorithm and binarized in the Hoechst channel, which was then used to generate a nuclear mask image. The mask images were then subjected to the Watershed algorithm in ImageJ to obtain single-cell regions of interest (ROIs). From the nuclei, the cytoplasm was captured by centripetal expansion of either of 3 pixels in images obtained with a 10X objective or of 6 pixels in images obtained with a 40X objective, until cell reaching the cell boundaries (cell membrane). The cytoplasm was then defined as the region between the cell membrane and the nucleus. Following cell segmentation, these cell boundaries were used to compute mean and integrated intensity values from all channels. Because ROIs are (initially) defined only by the nuclear signal, this approach is likely to over-or under-segment cells with irregular shapes, which can lead to nuclear, cytosolic or cell membrane 'signal contamination' between neighboring and/or stacked cells. Further experimental (e.g. including membrane markers to guide whole-cell rather than nuclear-only segmentation) and analytical algorithms to more accurately segment individual cells (e.g. using deep learning methods to register and apply additional features) would help to improve segmentation. All imageJ scripts used in this manuscript can be found in our Github repository (https://github.com/sorgerlab/cycif [Lin, 2018]; copy archived at https://github.com/elifesciencespublications/cycif).

Image stitching, shading and flat-field correlation
The BaSiC algorithm (Peng et al., 2017) was used for shade and flat-field correction in the create of the multi-panel montage images shown in Figures 2B, 6B, 9A and 11A. Additional information can be found on the BaSiC website (https://www.helmholtz-muenchen.de/icb/research/groups/quantitative-single-cell-dynamics/software/basic/index.html). An example of the performance of BaSiC is shown in Figure 2-figure supplement 1. The ImageJ plugin of BaSiC was applied for whole image stacks using the default options. After processing with BaSiC, images stack were stitched with ImageJ/Fiji 'Grid stitch' plugin with default options. ASHLAR was used to stich, register and scale images available at http://www.cycif.org/.

Time considerations
We believe that the greater time invested in t-CyCIF as compared to conventional IF IHC must be placed in the context of the much greater amount of data generate from a t-CyCIF experiment. It is also important to note that while t-CyCIF can be relatively slow when a single sample is processed it can easily be performed in parallel on multiple samples. As a practical example, we usually stain 30 slides in parallel (each involving 100-200 fields of view); in the case of TMAs, >80 samples can be assembled on each slide, so up to 2400 samples can be processed in parallel. With a single scanner, 30 slides can be scanned (average scan time~10 min) in about 6 hr. Photo-inactivation and washing steps take~1 to 1.5 hr, after which an additional round of staining is initiated. As a matter of convenience, we usually perform staining overnight. Hence, one user can generate data for 90 channels and 1800 images per day. Thus,~10 work days are required to generate 900 channels/18,000 images. Further time needs to be allotted for registration and stitching (~12-18 hr of computing time) and quantification (~24-48 hr computing time, depending on cell density). Overall, we believe that this is a reasonable level of throughput; moreover we have not yet attempted to optimize it using fluidic devices, automated stainers etc. We also note that the throughput of t-CyCIF compares favorably with other tissue-imaging platforms and single-cell transcriptome profiling.

Analysis of tissue integrity over cycles
We purchased a TMA (MTU481, Biomax Inc, https://www.biomax.us/tissue-arrays/Multiple_Organ/ MTU481) to test the impact of cycle number on tissue integrity. Images were captured and processed as described above. The registered image stacks were then segmented and nuclei counts for each core and each cycle were recorded. All values were normalized to the number of nuclei from the first cycle of a particular core biopsy and the fractional normalized nuclei count shown at each staining cycle.

Calculation of intensity overlap between different cycles and dynamic range
To compare staining patterns between different cycles within the same specimen, we calculated overlap integrals. First, we determined the distribution of intensity data averaged over each single cell and for each t-CyCIF cycles. The area under the curve of these distributions was calculated by trapezoidal numerical integration using 'trapz' function in Matlab (Gustafsson et al., 2008). The ratio of the area under the curve (AUC) for different cycles, samples or antibodies was calculated and the overlap scores then computed as: Overlap score ¼ overlap AUC=total AUC The dynamic range (DR) of fluorescence intensities for a given antibody was calculated as a rough estimate of the signal-to-noise ratio; SNR. The calculation was performed as follows: first, pixel-bypixel intensity data was extracted from a t-CyCIF image; the DR was then calculated as the ratio of the intensities of the 95th and 5th percentile values and represented on a log scale. High DR values indicate a favorable SNR. Intensities below the 5th percentile were considered to be background noise.

High-dimensional single-cell analysis by t-SNE
Raw intensity data generated from registered and segmented images were imported into Matlab and converted to comma separated value (csv) files. The viSNE implementation of t-SNE and EMGM algorithms from the CYT single-cell analysis package were obtained from the Pe'er laboratory at Columbia University (Amir et al., 2013). Intensity-based measurements (such as flow cytometry or imaging cytometry) of protein expression have approximately log-normal distribution (Bagwell, 2005), hence, t-CyCIF raw intensity values were first transformed in log or in inverse hyperbolic sine (asinh) using the default Matlab function or the CYT package (Amir et al., 2013), respectively. Between-sample variation was normalized on a per-channel basis by using the CYT package to align intensity measurements that encompass values between 1st and the 99th percentile. Data files were aggregated and used to generate viSNE plots. All viSNE/t-SNE analyses used the following settings: perplexity À30, epsilon = 500, lie factor = 4 for initial 100 iterations and lie factor À1 for remaining iterations.

Regional and neighboring analysis using K-nearest neighbors (KNN) methods
To determine whether PD-1 and PD-L1 expressing cells are sufficiently close for the receptor and ligand to interact, the spatial densities for PD1 + and PDL1 + cells were estimated using a k nearest neighbors (kNN) model with k = 4, corresponding to a~10 mm smoothing window. Since the density in space of the PD1 + or PDL1 + cells at any point in that space is proportional to the probability of that cell having a centroid there, the co-occurrence probability at a point was therefore proportional to the product of the spatial densities for both cell types at a point. To normalize for the difference in total PDL1+ or PD1+ cells between regions of the tissue corresponding to tumor and stroma, we calculated spatial probabilities for the different regions in the specimen separately. Figure 9-figure supplement 1 shows the distribution of co-occurrence densities for stroma and tumor relevant to a clear-cell carcinoma shown in Figure 9.

Calculating Shannon entropy values
Images were divided into regular grids and 1000 cells from each region used to calculate the nonparametric Shannon entropy as follows: Shanon Entropy s ð Þ ¼ Ài X s 2 i log s 2 i À Á where s i is the per-pixel intensity of signal s at a given point. Normalized Shannon entropy as calculated as E normalized = E region /E sample.

Expectation-Maximization Gaussian mixtures (EMGM) clustering
To determine an appropriate number of clusters (k) for analysis of the GBM tumor shown in Figures 11 and 12 and in Figure 12-figure supplement 2 we determined negative log-likelihoodratios for various values of k. For each choice of cluster number n, the likelihood-ratio was calculated for a Gaussian mixture model with n = k-1 and with n = k and the ratio then plotted relative to k. The EMGM algorithm was initialized 30 times for each value of k and it converged in all instances. The inflection at k = 8 (red arrow) suggests that inclusion of additional clusters (k > 8) explains a smaller, distinct source of variation in the data (Figure 12-figure supplement 1). As an alternative, k = 12 was also explored in Figure 12-figure supplement 2. Intensity values from all antibody channels (plus area and Hoechst intensity) were used for clustering.

Data availability
All data generated or analyzed during this study are included in the manuscript and supporting files. Intensity data used to generate figures is available in supplementary materials and can be downloaded from the HMS LINCS Center Publication Page (http://lincs.hms.harvard.edu/lin-elife-2018/) (RRID:SCR_016370).

Image availability
All images can be obtained from an OMERO image database via links found at the HMS LINCS Center Publication Page http://lincs.hms.harvard.edu/lin-elife-2018/ (RRID: SCR_016370). Stitched and registered image composites can be obtained at www.cycif.org. (RRID:SCR_016267) and via links found there. Data availability All data generated or analyzed during this study are included in the manuscript and supporting files. Intensity data used to generate figures is available in supplementary materials and can be downloaded from the HMS LINCS Center Publication Page (http://lincs.hms.harvard.edu/lin-elife-2018/) (RRID:SCR_016370). The images described are available at http://www.cycif.org/ (RRID:SCR_016267) and via and OMERO server as described at the LINCS Publication Page.