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

Abstract

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.

Keywords: cancer biology; computational biology; human; immunopathology; multiplexed imaging; single-cell method; systems biology.

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 Figure 1—figure supplement 1).

Figure 1—figure supplement 1.. Reduction in background signal intensity with repeated cycles of bleaching.

(A–C) Intensity distributions for three fluorescence channels (FITC/Alexa-488, Cy3/Alexa-555 and Cy5/Alexa-647) prior to pre-bleaching (blue), and after 1, 2 or 3 cycles of bleaching (black, red and green lines, respectively). With increasing number of bleaching cycles, the background signal is reduced 10- to 100-fold.

Figure 1—figure supplement 2.. t-CyCIF using antibodies labelled with Zenon Alexa-555 Fab fragments.

A human tonsil specimen was stained with unconjugated anti-CD11b antibody and then with Alexa-488 conjugated anti-Rabbit secondary antibody (left) or with the same anti-CD11b antibody following incubation with Zenon Alexa-555 (ThermoFischer; right). The Zenon Fab fragments generate non-covalent immune complexes that ‘label’ the primary antibody in a manner that is stable to subsequent processing steps.

Figure 2.. 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, α-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).

Figure 2—figure supplement 1.. Flat-field and shading correction for stitched images.

A large melanoma specimen (same as shown specimen as in Figure 2A–B) was imaged with a 20X/0.8NA objective on a CyteFinder, and a total 165 frame images were assembled to one image using the ASHLAR algorithm. Representative channels (Hoechst 33342 and S100-Alexa-488) of the stitched images before (left) and after (right) correction for uneven illumination using the BaSiC algorithm (Peng et al., 2017) (see Materials and methods).

Figure 2—figure supplement 2.. OMX super-resolution t-CyCIF images.

(A) Single-channel images for the composite image shown in Figure 2E. This melanoma sample was imaged using the GE OMX Blaze structured illumination microscope with 60X/1.42NA objectives. (B). Single-channel images for the composite image shown in Figure 2F. This breast cancer xenograft sample was also imaged using the OMX Blaze.

Figure 3.. t-CyCIF imaging of normal tissues.

(A) Selected images of a tonsil specimen subjected to 10-cycle t-CyCIF to demonstrate tissue, cellular, and subcellular localization of tissue and immune markers (see Supplementary file 1 for a list of antibodies). (B) Selected cycles from (A) demonstrating sub-nuclear features (Ki67 staining, cycle 1), immune cell distribution (cycle 2), structural proteins (E-Cadherin and Vimentin, cycle 5) and nuclear vs. cytosolic localization of transcription factors (NF-kB, cycle 6). (C) Five-cycle t-CyCIF of human skin to show the tight localization of some auto-fluorescence signals (Cycle 0), the elimination of these signals after pre-staining (Cycle 1), and the dispersal of rare cell types within a complex layered tissue (see Supplementary file 1 for a list of the antibodies).

Figure 4.. Efficacy of fluorophore inactivation and preservation of tissue integrity.

(A) Exemplary image of a human tonsil stained with PCNA-Alexa 488 that underwent 0, 15, 30 or 60 min of fluorophore inactivation. (B) Effect of bleaching duration on the distribution of anti-PCNA-Alexa 488 staining intensities for samples used in (A). The distribution is computed from mean values for the fluorescence intensities across all cells in the image that were successfully segmented. The gray band denotes the range of background florescence intensities (below 6.2 in log scale). (C) Effect of bleaching duration on mean intensity for nine antibodies conjugated to Alexa fluor 488, efluor 570 or Alexa fluor 647. Intensities were determined as in (B). The gray band denotes the range of background florescence intensities. (D) Impact of t-CyCIF cycle number on tissue integrity for four exemplary tissue cores. Nuclei present in the first cycle are labeled in red and those present after the 10th cycle are in green. The numbers at the bottom of the images represent nuclear counts in cycle 1 (red) and cycle 10 (green), respectively. (E) Impact of t-CyCIF cycle number on the integrity of a TMA containing 48 biopsies obtained from 16 different healthy and tumor tissues (see Materials and methods for TMA details) stained with 10 rounds of t-CyCIF. The number of nuclei remaining in each core was computed relative to the starting value; small fluctuations in cell count explain values > 1.0 and arise from errors in image segmentation. Data for six different breast cores is shown to the right. (F) Nuclear staining of a melanoma specimen subjected to 20 cycles of t-CyCIF emphasizes the preservation of tissue integrity (22 ± 4%). (G) Selected images of the specimen in (F) from cycles 0, 5, 15 and 20.

Figure 4—figure supplement 1.. Impact of bleaching time on fluorophore inactivation.

Intensity distributions for specimens stained with an antibody against Ki67 coupled to the following fluorophores: (A) Alexa-488, (B) Alexa-570 and (C) Alexa-647 prior to bleaching (blue curves) and after 15, 30 or 45 min of bleaching (black, red and green curves respectively). The distributions were calculated from the average fluorescence intensities of single cells as determined after image segmentation.

Figure 5.. Design of a 16-cyle experiment used to assess the reliability of t-CyCIF data.

(A) t-CyCIF experiment involving two immediately adjacent tissue slices cut from the same block of tonsil tissue (Slide A and Slide B). The antibodies used in each cycle are shown (antibodies are described in Supplementary file 2). Highlighted in blue are cycles in which the same antibodies were used on slides A and B at the same time to assess reproducibility. Highlighted in yellow are cycles in which antibodies targeting PCNA, Vimentin and Tubulin were used repeatedly on both slides A and B to assess repeatability. Blue arrows connecting Slides A and B show how antibodies were swapped among cycles. (B) Representative images of Slide A (top panels) and Slide B specimens (bottom panels) after each t-CyCIF cycle. The color coding highlighting specific cycles is the same as in A.

Figure 6.. Impact of cycle number on repeatability, reproducibility and strength of t-CyCIF immuno-staining.

(A) Plots on left: comparison of staining intensity for anti-PCNA Alexa 488 (top), anti-vimentin Alexa 555 (middle) and anti-tubulin Alexa 647 (bottom) in cycle 3 vs. 16 and cycle 7 vs. 12 of the 16-cycle t-CyCIF experiment show in Figure 5. Intensity values were integrated across whole cells and the comparison is made on a cell-by-cell basis. Spearman’s correlation coefficients are shown. Plots in middle: intensity distributions at cycles 3 (blue), 7 (yellow), 12 (red) and 16 (green); intensity 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 are available in a Jupyter notebook; see Code Availability section of Materials and methods for details.

Figure 6—figure supplement 1.. Comparison of staining intensities across different cycles at a single-cell level.

Images come from the antibody-swap and repeat experiment showing in Figures 5 and 6. Comparison of single-cell level intensities for (A) PCNA Alexa-488, (B) Vimentin Alexa-555 and (C) Tubulin Alexa-647 stained in four different cycles (cycle 3, 7, 12 and16). For any given cycle pair, single-cell intensities density plots are shown. Spearman’s correlation coefficients (rho) and the mean intensity ratios between two cycles are shown.

Figure 7.. t-CyCIF of a large resection specimen from a patient with pancreatic cancer.

(A) H&E staining of pancreatic ductal adenocarcinoma (PDAC) resection specimen that includes portions of cancer and non-malignant pancreatic tissue and small intestine. (B) The entire sample comprising 143 stitched 10X fields of view is shown. Fields that were used for downstream analysis are highlighted by yellow boxes. (C) A representative field of normal intestine across 8 t-CyCIF rounds; see Supplementary file 3 for a list of antibodies. (D) Segmentation data for four antibodies; the color indicates fluorescence intensity (blue = low, red = high). (E) Quantitative single-cell signal intensities of 24 proteins (rows) measured in ~4×10³ 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 Figure 8.

Figure 7—figure supplement 1.. t-CyCIF for examining large resection specimens of a human pancreatic cancer.

(A) Representative frame of small intestine from the PDAC resection shown in Figure 7 with images for PCNA, beta-catenin, Ki67 and pERk shown. These frames correspond to the segmented panels shown in 7D. (B) Creation of nuclear masks following identification of nuclei. Left panel: Hoechst image from t-CyCIF cycle one on the PDAC resection sample; middle panel: a binarized nuclear mask from cycle 8 (the final cycle in this t-CyCIF experiment); right panel: the overlay image of Hoechst stain (blue) and the cycle eight nuclear mask (yellow). The final cycle is used to create the mask used in this analysis so that the same cells can be tracked through all t-CyCIF cycles despite ~15% overall cell loss by cycle 8. Thus, regions in the overlay that show up in blue correspond in most cases to cells that are lost in the course of t-CyCIF and not to failure to identity and segment cells correctly.

Figure 8.. High-dimensional single-cell analysis of human pancreatic cancer sample with t-CyCIF.

(A) t-SNE plots of cells derived from small intestine (left) or the PDAC region (right) of the specimen shown in Figure 7 with the fluorescence intensities for markers of proliferation (PCNA and Ki67) and signaling (pERK and β-catenin) overlaid on the plots as heat maps. In both tissue types, there exists substantial heterogeneity: circled areas indicate the relationship between pERK and β-catenin levels in cells and represent positive (‘a’), negative (‘b’) or no association (‘c’) between these markers. (B) Representative frames of normal pancreas and pancreatic ductal adenocarcinoma from the 8-cycle t-CyCIF staining of the same resection specimen from Figure 7. (C) t-SNE representation and clustering of single cells from normal pancreatic tissue (red), small intestine (blue) and pancreatic cancer (green). Projected onto the origin of each cell in t-SNE space are intensity measures for selected markers demonstrating distinct staining patterns. (D) Fluorescence intensity distributions for selected markers in small intestine, pancreas and PDAC.

Figure 9.. Spatial distribution of immune infiltrates and checkpoint proteins.

(A) Low-magnification image of a clear cell renal cancer subjected to 12-cycle t-CyCIF (see Supplementary file 4 for a list of antibodies). Regions high in α-smooth muscle actin (α-SMA) correspond to stromal components of the tumor, those low in α-SMA represent regions enriched for malignant cells. (B) Representative images from selected t-CyCIF channels are shown. (C) Quantitative assessment of total lymphocytic cell infiltrates (CD3⁺ cells), CD8⁺ T lymphocytes, cells expressing PD-1 or its ligand PD-L1 or the VEGFR2 for the entire tumor or for α-SMA^high and α-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 µm of each other and with high spatial density (in yellow) and thus, are potentially positioned to interact at a molecular level.

Figure 9—figure supplement 1.. Spatial analysis of PD-1 and PD-L1 expressing cells.

This figure is relevant to the t-CyCIF study on a clear cell cancer shown in Figure 9. To compare co-localization of PD-1 and PD-L1 in tumor vs stroma, we computed the probability that PD1+ and PDL1+ cells would co-occur within a radius of ~10 um, normalizing for the difference in the total number of PD1+ and PDL1+ cells in the two tissue regions. We interpret the spatial density of PD1+ or PDL1+ cells at each point in space as proportional to the probability of their occurring there (see Figure 9E–F). The co-occurrence density at a point (Figure 9G) is therefore the product of the spatial densities for PD1+ or PDL1 +cells at that point. For simplicity, regions corresponding to tumor and stroma regions were defined by a diagonal line that separated the upper-right and lower-left regions of the tissue; this corresponded closely to α-SMA-low (tumor) and α-SMA-high (tumor) domains (left panel). We found an We found an e^(0.98)=2.7 fold difference in the distributions of co-occurrence densities between the two regions, representing an effect size of 0.94 as measured by Hedge's g; this represents a significant and potentially meaningful fold-change (right panel).

Figure 10.. Eight-cycle t-CyCIF of a tissue microarray (TMA) including 13 normal tissues and corresponding tumor types.

The TMA includes normal tissue types, and corresponding high- and low-grade tumors, for a total of 39 specimens (see Supplementary file 3 for antibodies and Supplementary file 5 for specifications of the TMA). (A) Selected images of different tissues illustrating the quality of t-CyCIF images (additional examples shown in Figure 9—figure supplement 1; full data available online at www.cycif.org). (B) t-SNE plot of single-cell intensities of all 39 cores; 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.

Figure 10—figure supplement 1.. Gallery of exemplary tissues imaged on the TMA described in Figure 10.

Gallery of exemplary tissues imaged on the TMA described in Figure 10.

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).

Figure 12.. Spatial distribution of molecular phenotypes in a single GBM.

(A) Clustering of intensity values for 30 antibodies in a 10-cycle t-CyCIF analysis integrated over each whole cell based on images shown in 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—figure supplement 2. (B) EMGM clusters (in color code) mapped back to the positions of individual cells in the tumor. The coordinate system is the same as in Figure 11A. The positions of seven macroscopic regions (R1-R7) representing distinct lobes of the tumor are also shown. (C) Magnified view of Frame 147 from region R5 with EMGM cluster assignment for each cell in the frame; dots represent the centroids of single cells. (D) The proportional representation of EMGM clusters in each tumor region as defined in panel (B).

Figure 12—figure supplement 1.. Determination of cluster number for semi-supervised clustering using expectation–maximization Gaussian mixture (EMGM) modeling.

To determine an appropriate number of clusters (k) for analysis of the GBM tumor shown in Figure 12 we determined negative log-likelihood-ratio for various values of k. Due to the large sample size, likelihood-ratio tests were not helpful in choosing k. Thus, 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) suggested that inclusion of additional clusters (k > 8) explains a smaller, distinct source of variation in the data. The plot is shown on a logarithmic scale to better visualize the range of the log-likelihood ratios, and should not be confused with the logarithm already applied to the likelihood ratios themselves.

Figure 12—figure supplement 2.. Spatial distribution of molecular phenotypes in a single GBM.

This analysis is directly analogous to the analysis shown in Figure 11, but uses 12 clusters for EMGM analysis rather than 8. (A) Intensity values from the tumor in Figure 10 were clustered using expected-maximization with Gaussian mixtures (EMGM) with k = 12. The number of cells in each cluster is shown as a percentage of all cells in the tumor. (B) EMGM clusters (in color code) mapped back to singles cells and their positions in the tumor. The coordinate system is the same as in Figure 10. The positions of seven macroscopic regions (R1-R7) representing distinct lobes of the tumour are shown. (C) Magnified view of Frame 147 from region R5 with EMGM cluster assignment for each cell in the frame shown as a dot. (D) The proportional representation of EMGM clusters in each tumor region as defined in Panel B.