The Normal and Fibrotic Mouse Lung in Situ Classied by High Dimensional Single Cell Analysis on Routinely Processed Tissue

Single cell classication is elucidating homeostasis and pathology in tissues and whole organs. We applied in situ spatial proteomics by multiplex antibody staining to routinely processed mouse lung, healthy and during a brosis model. With a limited validated antibody panel (24) we classify the normal constituents (alveolar type I and II, bronchial epithelia, endothelial, muscular, stromal and hematopoietic cells) and by quantitative measurements, we show the progress of lung brosis over a 4 weeks course, the changing landscape and the cell-specic quantitative variation of a multidrug transporter. An early decline in AT2 alveolar cells and a progressive increase in stromal cells seems at the core of the brotic process.


Introduction
Single cell biology, brought to fruition by advances in gene sequencing and computational progress, has revolutionized how we understand biological processes in health and in pathology 1 . .
Applying these techniques to the analysis of individual cells in-situ, i.e. within the tissue microenvironment, has added the information of the tissue sociology of the specimen, answering to the growing need to investigate it, due to various cellular functions, the spatial organization of molecular targets, the relationship among multiple cell types and morphology.
In situ high dimensional single cell analysis has been approached with genetic (in situ hybridization 2 ), proteomic (in situ MALDI-MS 3 ) and targeted phenotypic analysis via antibodies. To this latter end, the application of multiple antibodies to a single tissue section has been accomplished essentially by two methods: heavy metal isotope (IMC 4 ], MIBI 5 ) or uorochrome-tagged antibodies 6 . The uorochromes can be covalently linked to the antibodies or linked to DNA-barcodes and removed by sequential quenching (CyclIF 7 ), DNA lysis (CODEX 8 ), or physical removal together with the antibody from the section 9 .
While in situ single cell sequencing can provide information on the RNA species in the thousands, in situ multiplex staining has a limit of around 100 antibodies, in general 30-60 for the time being 6 . To distinguish the techniques which can add antibodies in the dozens from the other which provide a much limited number (5-6), we will refer to the former as "high-plex" multiplexing techniques.
Single cell biology via high-plex in situ staining can be accomplished both in fresh frozen tissue sections and in sections from routinely processed tissue (formalin xed, para n embedded; FFPE), depending on the technology used. FFPE material has several advantages, including the capability to create high density arrays containing several different samples, the tissue microarray technology (TMA) 10 . Being the high-plex technology complex and time-consuming by itself, the combination with the TMA technology allows an enormous, detailed analytical power in the basic and applied science eld.
Murine models are an extremely helpful, standard companion to the research on human subjects in all elds of medicine. All kinds of diagnostic and therapeutic approaches used on humans have been miniaturized and applied to mice 11 , including histopathologic and phenotypic examination. The problem of using antibody-based high-plex methods on mouse tissue when using primary antibodies raised in mice is the cross reactivity with endogenous immunoglobulins: we have overcome this problem by using anti-isotype secondary Abs 12 on FFPE mouse tissue, therefore we have broadened the antibody portfolio which can be used on mouse histopathology.
Mouse specimens have been extensively used for single cell biology studies in all elds, largely via single cell RNA sequencing, rarely for in situ spatial proteomics and high-plex staining.
We have applied the MILAN high-plex technique 13 to a murine model of bleomycin-induced (BLM) lung brosis. This animal model is the most used in preclinical studies 14 and, in this context, the American Thoracic Society (ATS) has suggested guidelines to ensure that in vivo animal modeling studies have the highest chance of discriminating between potentially effective and ineffective anti brotic compounds [15][16][17] .
Bleomycin elicits a time-dependent lung brotic process 14 , which has several similarities with the human counterpart 18 . To better understand lung brosis models, highplex technology has been used to detect any target of interest and apply to routinely processed samples with minimal requirements.
In this study we consider the recent observation that BLM administration induces an ABC (ATP Binding Cassette) carrier upregulation in lung, with an increasing expression of P-glycoprotein (P-gp) and transporters in C57BL/6 male mice (Park et al 2020) 19 . The ABC carriers actively transport multiple xenobiotics across the membrane reducing the intracellular concentration of drugs and leading to a potential decrease in anti-brotic activity.
We present here the validation of the method and results of the single cell composition of the healthy and pathologic mouse lung, including the changes in expression of a P-gP multidrug transporter in each of the cell types identi ed by high-dimensional analysis.

TMA validation
The use of tissue cores, instead of whole sections, entails higher throughput at the expense of representativeness. To address this latter aspect we used single cell analysis to validate the use of the TMA. Single cell analysis included an average of 2396 cells for 1 mm core, 10337 for 2 mm core and 40904 for the whole section. By high dimensional analysis, the number of clusters found were respectively 12, 12, 13, 13 (1mm core), 17, 15, 17, 18 (2 mm core) and 22 (whole slide). Clusters comprising sparse cells, such as Macrophages (Fig. 1) and T cells (not shown) were equally represented in the 2mm cores and in the whole section (Fig. 1).
The 1 mm cores, instead, poorly identi ed speci c cell populations (macrophages, T lymphocytes, etc.), thus it was not representative; 2mm TMA cores were used for the study. SINGLE CELL ANALYSIS OF THE NORMAL AND FIBROTIC MOUSE LUNG tSNE plots from four control animals showed a superimposition of all cells from all animals as expected (Fig 2A). Cells from bleomycin-treated mice were allocated in the same phenoclusters (Fig 2E) of the controls, but with a slightly different spatial distribution for mice 7, 14 and 21 post BLM treatment ( Fig  2B). Mice at day 28th were very similar to controls (Fig 2C). The phenoclusters contained cells from every case (Fig. 2D) demonstrating the absence of batch effect in the data.
Subsequently, all cases were analyzed and each cluster classi ed individually ( Supplementary Fig. 4) After the removal of clusters resulting from artifacts or uninterpretable clusters (see supplementary material), a total of 80506 cells were analyzed in the 13 cores considered, with an average distribution of 6029 cells/core. The normal mouse lung contained a majority of Alveolar Epithelial cells (67.9% ± 11.9), divided into AT1 (22% ± 5.8), AT2 (14% ± 9.7) and a population with coexpression of AT1 and AT2 markers, named transitional AT (32% ± 11.8). Each of the remaining constituents of the normal lung remained below the 10% value. Bronchial epithelial cells and goblet cells averaged 7% ± 7.7, endothelial cells were 3% ± 1.4, smooth muscle cells 2% ± 1.9, stromal cells 0.7% ± 0.4. Resident hematopoietic cells were T lymphocytes (7% ± 1.5), B lymphocytes (3% ± 1.9), macrophages and dendritic cells (4% ± 1.9) and neutrophils (1% ± Single cells from mice treated with Bleomicyn and examined at day 7, 14, 21 and 28 partially overlap with the control mice (Fig 2), highlighting qualitative and quantitative changes in the lung population. By applying the same classi cation criteria used to classify the normal lung, one could observe a transient reduction in AT1 and a progressive decrease of AT2 pneumocytes (Fig. 3), largely due to changes in an AT2 subpopulation (transitional AT2, see supplementary Fig. 5 and 6). Stromal cells increased progressively, accompanied by a late increase of smooth muscle cells, mimicking the histopathologic accumulation of stroma (see supplementary Fig. 1). The remaining lung cells remained stable except for a late decrease of macrophages, neutrophils and B lymphocytes. Plasma cell markers were not included in the panel, thus we could not assess whether the decrease in B cells was due to absence or maturation.
The spatial distribution of the cell types revealed a crowded parenchyma and focally increased stromal foci, B cell and macrophage aggregations (Fig. 4).
A drug transmembrane transporter, P-gP, constitutively expressed in many cell types, was measured in lung cells identi ed by high dimensional analysis. Increasing intensity signal was registered in all populations after bleomycin treatment, except in bronchial cells. The onset of the time-dependent increase was delayed in the stromal cell type (Fig. 5);

Discussion
We have obtained an in situ spatial representation of the normal and of the brotic mouse lung by single cell classi cation via multiple antibody staining. Analogous to a single cell classi cation of murine (and human) lung cell component by single cell RNA sequencing 20-23 , here we show that the content of a normal or diseased organ such as the lung can be nely dissected at the single cell level, with two additional properties: spatial cell disposition is represented and proteins, instead of RNA, are assessed by a robust, cheap and versatile method.
With a rather limited number of validated antibodies, we can identify in situ the main cell types which are relevant for lung homeostasis and for the initiation and establishment of the bleomycin-induced brosis: AT1, AT2, transitional AT1-AT2 pneumocytes, bronchial lining cells, vasculature, stromal cells, T and B lymphocytes, macrophages. In addition we measured at the single cell level continuously expressed proteins such as a transcellular drug transporter, P-gP.
We measured with the same detail the changes occurring upon Bleomycin treatment over time, documenting a change in the AT1/AT2 ratio, the progressive accumulation of stromal cells and the asynchronous increase of P-gP expression.
Lastly, we included in the assay a high throughput method, the tissue microarray technology, after careful validation of the representativeness of tissue core size.
We suggest the use of a 2 mm TMA core as a minimum tissue area for the evaluation of tissue cell composition on mouse lung. We came to this conclusion by performing high-dimensional analysis and cell clustering as a validating tool, applicable to TMAs from tissues of various origin and composition.
The most abundant population we found in normal tissue and also in BLM-treated mice is represented by epithelial cells.
Data which has been obtained by tissue dissociation and single cell sequencing of normal mouse lungs 24,25 do not agree among themselves about the representation of the cell types, probably because of different pre-analytical dissociation methods and inherent selective loss of epithelial cells or enrichment of other cells. The drawbacks of the tissue dissociation methods are known 26 . By in-situ cell classi cation we can provide an unbiased estimate of the various lung cell components. In addition, by measuring the end product of the transcription and translation machinery, we complement the data provided by RNA sequencing.
Analogously to scRNA sequencing data, the alveolar epithelial cell populations do aggregate in closed contact in the bidimensional tSNE space, re ecting a continuum of phenotypes, interpretable as a transition from AT2 to AT1 cells. We tentatively identi ed a reproducible subset of AT2 cells bearing AT1 markers, which we dubbed "transitional AT" cells, analogously to a similar subset identi ed in Bleomycintreated animals 20,22 . Whether the subset we identi ed by in-situ proteomics and the ones described in scRNAseq experiments are identical requires additional experiments with parallel analytical tools.
Collagen deposing stromal cells are the main actor in the brotic process produced in this experimental model. We could demonstrate an increase in stromal cells which parallels the histomorphologic changes. According to our data, these cells represent a minority of the lung population, to the point that are inconsistently demonstrated in the normal lung, at variance with data obtained by dissociation and RNA sequencing. Besides the differences in methods, as outlined before, stromal cells are underrepresented in our markers panel; in addition, a nuclear-based cell segmentation is not ideal to identify elongated cells, thus we may underestimate this cell subset, despite showing a treatment-dependent increase.
Histochemical stains do not discriminate the cell of origin for collagen deposition, nor if the collagen is deposited by pre-existing cells. On the other hand we assess stromal cells individually and by intracellular markers, documenting a net increase in nucleated stromal cells upon treatment: it is thus not unexpected that the two represent non-identical assays, as demonstrated by the data.
The progressive decrease of the AT2 and particular of the transitional AT over time and the increase of stromal cells we have shown may highlight the key drivers of lung brosis in this mouse model: a progressive reduction of the alveolar lining repair by depleting local alveolar progenitors and surfactant producing cells, AT2, coupled with newly produced, collagen depositing broblasts. These results are novel and should be consolidated by independent experiments and increased sample numerosity.
The time-course expression of P-gP has been previously published for the same Bleomycin-induced mouse lung model 19 , with notable differences: the experiments were performed on male mice and the RNAs for the multidrug transporters were measured, both as whole tissue extracts and by in-situ hybridization. We have reproduced the data on male mice (see Supplementary data), and shown that female mice differ in the kinetic of the pathology but not in the type of histopathologic lung changes. In addition, by measuring the multidrug transporter protein at the single cell level on each and every cell type in the specimen, we have detailed a cell-type speci c kinetics of P-gP expression, with some notable difference with the RNA data of Park et. al. 19 .
Limitations of the study: We employed a limited number of antibodies, compared to the potentiality of the MILAN technique (over 100); this because there is a variety of reagents to choose from which is more limited than what is available for human FFPE tissue, particularly in terms of species origins of the antibodies. Being the MILAN technique based on multiple (3 or 4) unconjugated antibodies per round, we partially overcome these limitations by using mouse antibodies on mouse tissue, background free 12 . It has been shown previously however, that the discriminating power of the dimensionality reduction algorithms is so high, that even with a reduced number of diagnostic parameters, the cell types of interests can still be identi ed 27 and this is the case for the present work.
An additional limitation of this study is the use of a rather unsophisticated image segmentation technique and manual cell type assignment for the phenogroups. Cell segmentation is one of the most challenging tasks for in situ transcriptomics and proteomics 28 and efforts are ongoing to improve it for the mouse lung brosis model. In this model, elongated cells such as broblasts and stromal cells in general poorly t cell identi cation via nuclear DAPI identi cation. This may be the reason why in normal lungs we occasionally found minimal or no broblast cell clusters, together with the small total amount, in common with cell suspension based studies [20][21][22][23] . Along the same lines, we found clusters containing markers of two distinct cells (e.g. epithelial and macrophages, epithelial and stroma) which we were forced to discard, because of the juxtaposition of elongated cells was not solvable with the cell segmentation strategy used, despite being clearly identi able by the tissue spatial distribution (see supplementary data/discussion).
In summary, we have shown proof of principle that mouse tissue, being a normal organ, a pathology model or a developing tissue, can be dissected on routinely processed material by in-situ highdimensional proteomics and single cell bioinformatic analysis. This represents a powerful tool for preclinical studies such as drug discovery and novel treatments and can be integrated with other "omics" tools such as scRNA sequencing, in situ transcriptomics etc. 29 .

Experimental animals and histology
A mouse model of bleomycin-induced lung brosis, currently in use in our institution, has been previously published 30  In brief it consists of the oropharyngeal BLM administration to 7-8 weeks old C57Bl/6 female mice of 15 µg/mouse at each day of treatment respectively. Animals were sacri ced at 7, 14, 21 and 28 days after the administration.
After the sacri ce, lungs were removed, formalin xed and para n embedded (FFPE). Three serial sections, 5 μm thick were obtained for the stainings: Hematoxylin and Eosin, Masson's trichrome (TM) and a multiplex immunostaining.
Whole-slide images were acquired by a NanoZoomer S-60 Digital slide scanner (NanoZoomer S60, Hamamatsu, Japan) at 20x magni cation. Fibrotic lung injury was assessed histologically through the Ashcroft scoring system 31,32 on whole parenchyma by two independent researchers (blinded to the experimental design). Detailed description and histomorphometric characterization are included in the Supplementary methods and in Supplementary Figure 1. The most relevant areas on the slide were marked for subsequent tissue microarray construction (see below).

Tissue Microarray design and construction
For the TMA preparation, we used FFPE tissue blocks of lungs from BLM-treated female mice. 13 cases of treated and control lungs were selected; cores were extracted, selecting speci c regions representing normal lung tissue with alveolar parenchyma, bronchioles and vessel for saline samples and broproliferative foci with its alterations for bleomycin ones from the C57BL/6 female mice FFPE lungs of BLM and saline treated mice across four time points (7,14,21 and 28 days, see Supplementary Table  1). The TMA was constructed with a Tissue Microarrayer Galileo CK4500 (Tissue Microarrayer Model TMA Galileo CK4500; Integrated Systems Engineering srl, Milano, Italy) using Galileo Software to match the annotated tissue on histological slide with their corresponding areas on the surface of the para n donor blocks. After the core transfer in the recipient block, to allow the samples to be properly embedded into the block, the TMA was incubated for 24h at +38°C. Finally, 5 μm thick serial sections were cut from the TMA with a rotary microtome (Slee Cut 6062, Slee Medical, Mainz, Germany) and placed on Polysine adhesion glass slides (Thermo Fisher Scienti c). Routine histology were performed as described previously 30 .
In order to select the appropriate Tissue Microarray (TMA) core dimension, after applying a high-plex MILAN staining method 13 on the whole slide of a mouse lung, we selected 4 virtual cores (with ImageJ) of 1 mm and 4 of 2 mm in diameter. Thus we run single cell analysis as previously described on 8 cores and on a whole slide and we compare the results. The level of adequacy of tissue portion was set when all main cell populations found in whole tissue by clustering analysis (Rphenograph package ) were rediscovered in cores.

Highplex Immuno uorescence
Lung samples were processed as per the Milan (Multiple Iterative Labeling by Antibody Neodeposition) protocol 13 modi ed for mouse FFPE tissue staining. The protocol consists in the cyclic application of primary antibodies (Table 1) raised in multiple species, including mouse 12 , applied on the same section, nuclear staining with DAPI, auto uorescence subtraction 9 . For each round, the stained slides were scanned and saved as digital images using a multichannel uorescence acquisition instrument (NanoZoomer S60, Hamamatsu, Japan). The antibody stripping preceded a subsequent staining round. A multichannel uorescence acquisition was performed after each stripping, in order to check the complete antibody removal.
In order to generate an antibody panel which would produce single cell lung classi cation, we mined existing single cell RNA sequencing (scRNAseq) studies 24,25 for highly expressed, lineage-restricted messages which would correspond to proteins against which antibodies would be available. The staining sequence design and antibodies speci cations were provided in Table 1. The methods used for the antibody selection and validation are reported in Supplementary methods.
Sections were incubated overnight with primary antibodies, which were then revealed by secondary uorochrome-tagged antibodies; both isotype control and omission of primary have been performed as negative control. Sequential stripping was obtained with a beta-mercaptoethanol and sodium dodecyl sulphate mix 9 .
Mouse tissue preservation after each of multiple staining and stripping rounds showed no or negligible cell loss ( Supplementary Fig 2 and 3), as measured by DAPI nuclear stain and quanti cation on four other tissues (kidney, liver, lung and heart).

Single cell analysis
After the stainings were acquired, digital slide images (.ndpi) were imported as .tiff and registered with the AMICO software 33 , based on Fiji. All DAPI images from different staining rounds were registered together and their coordinates of rotation and translation were used to align the individual marker images of each round. Once all images were aligned, auto uorescence was subtracted from FITC and TRITC channels.
DAPI stainings were used for cell segmentation with Fiji (Threshold Default -Watershed_Analize particles and Count mask creation) Mean intensity value of all markers within each segmented cell and spatial coordinates of centroids of nuclei were recorded together in a .csv le.
Then .csv les were uploaded to the R Studio software (version 1.4) for a more detailed analysis. Rtsne 34 , umap 35 and Rphenograph 36 (n=30) algorithms (R packages) were used respectively for dimensionality reduction and clustering of data. UMAP plots were decorated with single individual relevant marker intensity to evaluate the distribution among phenogroups.
A single tSNE plot from homogeneous groups of animals was used to compare treatments vs control, then the samples were analyzed individually.
Clusters obtained were further explored via: i) individual cell immunopro ling through dedicated hierarchical clustering heatmaps as published 37 , ii) visualization of the spatial distribution on tissue. Clusters which satisfy both criteria were manually classi ed according to the expressions of key markers (Table 1).
Major cell types were created by merging clusters with a similar phenotype, after removing artifacts. Main types were organized as follows: 1) Alveolar Epithelial cells subdivided in AT1, AT2 or transitional AT1-AT2 subpopulations, 2) Bronchial and Goblet Epithelial cells, 3) Macrophage cells, 4) B cell, 5) T cells, 6) Vasculature and 7) Stromal cells, Neutrophils (8). For graphic spatial representation, all epithelial cells were grouped in a single entity (Alveolar Epithelial Cells).
Phenoclusters for which no coherent phenotype was identi able were excluded from the analysis (see supplementary material for further information).
The transporter expression was analyzed by measuring the mean of the intensity of the signal (on 8 bit grayscale images, values 0-255) within each main type. Tables   Table 1: Design and antibodies speci cations. Figure 1 See image above for gure legend