Samuel Kubler

Suvadip Mukherjee

Jean-Christophe Olivo-Marin

Thibault Lagache

A ROBUST AND VERSATILE FRAMEWORK TO COMPARE SPIKE DETECTION METHODS IN CALCIUM IMAGING OF NEURONAL ACTIVITY

Keywords: Calcium imaging, neuronal activity, simulation, deconvolution, Hydra

Calcium fluorescence imaging enables real-time activity monitoring of single neurons in living animals. A critical inverse problem resides in the precise inference of spike locations from noisy fluorescence traces, especially in the presence of burst spiking and non-linear fluorescence intensity. Several spike extraction algorithms have been proposed in the recent years, but a robust and objective evaluation of their performance still remains elusive due to the unsupervised nature of the problem. Here we propose a biologically-inspired mathematical framework to reproduce synthetic fluorescence traces from a time-series data of spike-trains. The idea is to create a versatile platform to objectively test the state-ofthe art spike inference methodologies over a large range of experimental parameters. Our solution appears as a complementary but more exhaustive approach to determine the robustness of existing solutions to different nature of signals, imaging artefacts, sensitivity to hyper-parameters and pre-processing steps. We benchmark state-of-the-art algorithms with the proposed simulation platform, and validate the results on an experimental dataset of the Hydra Vulgaris. We hypothesize that, in contrast to the common practice of qualitative evaluation, quantitative measure of algorithm robustness is essential in understanding the suitability of a spike inference algorithm to be used in an automated computational pipeline to decipher the neural code.

INTRODUCTION

To understand the emergent computational properties of connected single neurons, it is necessary to monitor the coordinated activity of many single cells with high temporal resolution [START_REF] Buzsáki | Large-scale recording of neuronal ensembles[END_REF]. Fluorescence calcium imaging [START_REF] Yuste | Control of postsynaptic ca2+ influx in developing neocortex by excitatory and inhibitory neurotransmitters[END_REF] remains the gold-standard for such studies. Thanks to the recent advances in calcium probe engineering and optical microscopy, thousands of interconnected neurons can now be imaged in living animals with high temporal resolution [START_REF] Yang | In vivo imaging of neural activity[END_REF]. The success of reverse engineering the brain now hinges upon effective signal processing techniques to extract critical information from * samuel.kubler@pasteur.fr imaging datasets. In addition to the automatic tracking of single neuron fluorescence [START_REF] Thibault Lagache | Robust single neuron tracking of calcium imaging in behaving hydra[END_REF], another critical issue in the processing of calcium imaging data is that of analysing neuronal spike trains from the videos of firing neurons (see Fig. 1). The inverse problem of reconstructing the neuronal spikes from the calcium imaging data is inherently unsupervised. Given a temporal sequence of calcium uptake of the firing neurons, it is difficult, even for a neuroscientist, to reconstruct a precise representation of the spiking pattern. State-of-the-art deconvolution methodologies exhibit significant variability in extracting the spike patterns from same videos [START_REF] Deneux | Accurate spike estimation from noisy calcium signals for ultrafast threedimensional imaging of large neuronal populations in vivo[END_REF] [6] [7]. Due to the absence of ground-truth annotation in most experimental data, qualitative evaluation appears to be a complementary norm to measure the efficacy of deconvolution algorithms [6]. However, supplemented by quantitative experimental evaluations, such qualitative assessments are unreliable due to inherent subjective bias. In the absence of gold standards, the proper mechanism to study the behaviour of deconvolution algorithms is via rigorous testing on realistic synthetic data.

In this paper we present a robust generative model to realize time-series distributions of fluorescent calcium traces from a given neuronal spike-train. Existing fluorescent trace simulators [START_REF] Deneux | Accurate spike estimation from noisy calcium signals for ultrafast threedimensional imaging of large neuronal populations in vivo[END_REF][START_REF] Giovannucci | Caiman: An open source tool for scalable calcium imaging data analysis[END_REF] do not always account for the intricate features which exist in real calcium imaging datasets. In contrast, our proposed methodology enables a realistic simulation by integrating four critical design factors: (a) a model for photobleaching of calcium indicators, (b) imaging noise model, Existing deconvolution techniques may be categorized in three major groups (see table 1): 1) deterministic methods that compute the optimal spike time series together with the estimated calcium dynamic of the fluorescent indicator [START_REF] Pnevmatikakis | Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data[END_REF];

2) probabilistic approaches that compute the convolved time series of spikes that maximises the likelihood of observed data [START_REF] Deneux | Accurate spike estimation from noisy calcium signals for ultrafast threedimensional imaging of large neuronal populations in vivo[END_REF]; and 3) machine learning approaches [6]. We study the robustness of four representative methodologies and evaluate to which extent they can be adapted on heterogeneous simulated signals with strong underlying dynamic. Finally, we show that the results provided by the simulator, and those obtained within experimental dataset from Hydra Vulgaris are correlated.

METHOD

Our simulator can be decomposed into four main parts, shown in Fig. 2 : 1) the modeling of the spike impulse train using the Firing Instantaneous Rate (FIR) function; 2) the modeling of the calcium trace from dirac impulse signal using a kernel learned from experimental dataset; 3) the modeling of the fluorescence trace from the calcium by handling the photobleaching effect through a Poisson-Gaussian noise; 4) an additive non-linear baseline integration. Each step is detailed hereafter.

To model the firing pattern s(t) of a single neuron, we implement a Poisson process with a FIR describing the probability at each time step that the neuron fires. The FIR is modeled as a piece-wise step function whose high constant steps correspond to burst events, and low steps to the background spontaneous firing activity.

FIR(x) = n i=1 γ i χ [ti,ti+di] (x). ( 1
)
Here χ is an indicator function, equals to 1 when x ∈ [t i , t i + d i ] and 0 otherwise. The firing-rate γ i follows two regimes (bursting event or not), and is modeled with a truncated Gaussian law to keep frequencies positive. The number of steps n is linked to the burst-rate parameter. The burst duration d i follows a Poisson law, and the spike burst locations t i are drawn uniformly on the available time interval. The spike impulse signal s(t) is then derived from the FIR using an adaptation of inhomogeneous Poisson process simulation by thinning [START_REF] Lewis | Simulation of nonhomogeneous poisson processes by thinning[END_REF]. An exponential course of return to the equilibrium is used after each firing (the observed collapses of the FIR shown in Fig. 2) to model calcium reservoir depletion.

c(t) = s(t) * k(t) with k(t) = Ae -( t τ D ) β 1 -e -t-µ τ R (2)
Here τ D is the time constant of the calcium concentration return to steady-state, β > 0 is a power law, µ is the median time of calcium increase after the electrical spike. τ R is the corresponding rising time constant. In our simulations, we calibrate the convolution kernel with data extracted from Hydra calcium traces. In the third step, we model the photo-bleaching of calcium indicators with a mono-exponential decrease [START_REF] Jezierska | Poisson-Gaussian noise parameter estimation in fluorescence microscopy imaging[END_REF] λ(t) = c(t)e -t τ , with τ being the photo-bleaching time constant. To account for the Poisson shot noise of microscopes, we further model the recorded signal as a Poisson process such that P (t) ∼ P oisson(αλ(t)), where α is the gain of the microscope [START_REF] Chenouard | Avancées en suivi probabiliste de particules pour l'imagerie biologique[END_REF].

To obtain a realistic calcium trace, we add a Gaussian noise G ∼ N (m, σ 2 ) with constant mean m and standard deviation σ to the recorded Poisson signal (mixed Poisson-Gaussian representation) and we also add a periodic deterministic baseline B(t) = Asin(2πf t) where A and f are tunable amplitude and frequency respectively. Finally, the fluorescence calcium trace y(t) is derived from:

y(t) = P (t) + G + B(t) (3)
The photobleaching time constant and additive Gaussian noise parameters are fitted to experimental datasets using least square method. In our simulations, the SNR of generated calcium traces is modulated using the gain α of the microscope detector [START_REF] Chenouard | Avancées en suivi probabiliste de particules pour l'imagerie biologique[END_REF], an increased gain leading to higher SNR.

RESULTS

Benchmarking Deconvolution Algorithms

Using synthetic fluorescence traces, we benchmark four despiking methods (two of them are issued from the CaImAn library [START_REF] Giovannucci | Caiman: An open source tool for scalable calcium imaging data analysis[END_REF]): 1) Deterministic OASIS [START_REF] Friedrich | Fast online deconvolution of calcium imaging data[END_REF]; 2) Deterministic CDfoopsi (Constrained Foopsi) [START_REF] Pnevmatikakis | Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data[END_REF]; 3) Probabilistic MLspike [START_REF] Deneux | Accurate spike estimation from noisy calcium signals for ultrafast threedimensional imaging of large neuronal populations in vivo[END_REF] and 4) a Naive method that consists of smoothing the signal with a wavelet thresholding, before computing the first derivative of the signal and estimates the spikes locations with derivatives greater that one standard deviation of the derivative over the entire calcium trace. After having benchmarked the methods using different pre-processings proposed in the literature such as low-frequency filtering and normalization [START_REF] Vogelstein | Fast Nonnegative Deconvolution for Spike Train Inference From Population Calcium Imaging[END_REF] or polynomial detrending [START_REF] Nolan | Algorithms for the correction of photobleaching[END_REF], we conclude that the efficiency of a pre-processing method depends on each deconvolution algorithm, and we choose to perform a polynomial detrending and data normalization for OASIS, and no preprocessing for CDfoopsi and Mlspike.

Comparing the accuracy and robustness of methods

Deterministic methods (OASIS and CDfoopsi) do not handle non-linear baseline. Therefore, high baseline amplitude A will lead to false spike burst detections. To tackle this technical issue, we filtered the inferred spikes with respect to their estimated amplitude using a user-defined Decisional Amplitude Threshold (DAT). The DAT is the only parameter we optimize in tested deterministic methods since the others have been exhaustively analysed in previous study [7]. For Naive thresholding method, the only parameter is the threshold (typically a multiplicative factor of the derivative standard deviation). The probabilistic method MLspike that concomitantly estimate the spike locations and the underlying fluorescence baseline presents four hyper-parameters: the relative amplitude a of the spikes, the decay time constant τ of the calcium fluorescence trace after a spike, the a priori level of Gaussian noise σ and the drift d of the estimated stochastic baseline. To compare the performances of the different deconvolution methods over a large range of simulation parameters, we compute for each method the rate of missing and false positive spike detections and summarize detector performances using the Error Rate (ER) indicator proposed in [START_REF] Deneux | Accurate spike estimation from noisy calcium signals for ultrafast threedimensional imaging of large neuronal populations in vivo[END_REF]. We first evaluate the accuracy of each method by varying the different parameters of our synthetic simulator (see Fig. 3) : the gain α that modulates the SNR, the rate of burst events, the spiking rate inside bursts and the baseline amplitude A and frequency f (n = 10 simulations per set of simulation parameters). The method hyperparameters are calibrated using a grid-search which minimizes the error rate with the available ground-truth. Each method calibration is performed once for a reference set of simulator parameters: Gain of the detector α = 1, Firerate In Burst f IB = 0.75 Hz, Firerate Out of Burst f OB = 0.15 Hz (firerates correspond to γ i in/out burst regime respectively see Eq. 1), Burstrate B = 0.02 Hz, Baseline Amplitude A LF = 50, Baseline frequency f LF = 2e -3 Hz (see Eq. 3). These values have been chosen to make the simulator relevant with the experimental Hydra Vulgaris dataset. We find that, after hyperparameter calibration, the probabilistic MLspike is overall the most accurate method for each set of parameters. However, we observe that MLSpike performance rapidly degrades as simulation parameters change, especially for the gain α (SNR) and baseline amplitude. The other methods are overall much more robust to parameter variations while being less accurate for the specific set of calibration parameters. We highlight that the naive derivative thresholding produces good results for isolated spikes but can not properly handle bursts. CDfoopsi method [START_REF] Pnevmatikakis | Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data[END_REF] seems to provide the best compromise between accuracy and robustness to parameter variations.

Correlating simulation results with Hydra experimental dataset

The trade-off between accuracy and robustness of a despiking method is important for experimental in vivo applications as different neurons can present different firing patterns (isolated or burst) and SNR for example. We therefore compare despiking methods using experimental dataset composed by ∼ 250 single neurons imaged over 2000 frames at 10 Hz inside freely-behaving animal Hydra [START_REF] Dupre | Non-overlapping Neural Networks in Hydra vulgaris[END_REF]. This dataset was obtained thanks to the robust monitoring of single neuron activity with tracking algorithm EM C 2 [START_REF] Thibault Lagache | Robust single neuron tracking of calcium imaging in behaving hydra[END_REF]. Method parameters (DAT for deterministic methods and hyper-parameters for MLspike) were calibrated by manually labelling spike locations on real fluorescence traces, and by minimizing the distance between the inferred spikes and the manual annotation. Hydra calcium traces are heterogeneous in terms of baseline, noise and spiking dynamics, with neurons spiking quite sparsely and other in bursts (Fig. 4, [START_REF] Dupre | Non-overlapping Neural Networks in Hydra vulgaris[END_REF]). As observed with previous simulations, the second deterministic method [START_REF] Pnevmatikakis | Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data[END_REF] is the most robust to calcium trace heterogeneity and presents the best average Error Rate compared to manual labelling (OASIS: ER avg = 0.40, CDfoopsi : ER avg = 0.18, Naive : ER avg = 0.28, MLspike : ER avg = 0.42) The average poor performance of MLspike is likely due to its poor robustness to changes in experimental conditions (see MLspike in Fig. 5). The lack of ground truth for validating the deconvolution methods on experimental dataset imposes to analyse the consistency of the inferred spike locations. This evaluation hinges on biological assumption of existing neuronal ensembles (e.g. RP and CB neural networks in Hydra [START_REF] Dupre | Non-overlapping Neural Networks in Hydra vulgaris[END_REF]) that are supposed to provide correlated spiking patterns. Also, a sufficient and balanced average number of spikes per trace should reflect the robustness of each method to heterogeneous baseline variation and noise level. In this regard, the spike distributions over the dataset are summarized in the raster-plots (see Fig. 5-a). We also calculate the average number of spikes per neurons (OASIS : 12.61 ± 4.43, CDfoopsi : 22.62 ± 14.70, Naive : 19.22 ± 4.47, MLspike : 17.51 ± 26.10) and the neuron pair-wise correlation matrices using a Jaccard distance relaxed by a time tolerance (see Fig. s(n * i ) is the binarized impulse spike signal inferred for neuron n * i . The heterogeneity in spike distribution varies between the methods. MLSpike computes many false spike detections (horizontal lines in raster-plot combined with a high standard deviation of the average number of spikes per sequence), while the naive thresholding method misses more spikes, but provides more balanced impulse spike signals in terms of average number of spikes per fluorescence trace. We observe that the highest global average Jaccard distance is obtained for the naive thresholding and CDfoopsi methods (see Fig. 5b). Finally, we conclude that CDfoopsi appears as the best tradeoff solution since it provides the best reconstruction results compared to our manual labelling, a sufficient and balanced number of spikes per trace and one of the highest global Jaccard correlation. The results obtained on experimental Hydra Vulgaris dataset are congruent with the ones derived from the proposed simulator which enforces the conclusion drawn by [7] about the versatility of CDfoopsi, especially within dataset with heterogenous baselines.

CONCLUSION

In this paper, we provided a mathematical framework for a fluorescence trace simulator to objectively compare and validate spike deconvolution techniques. State-of-the-art deconvolution methods were benchmarked on synthetic and experimental datasets. Such comparative analysis is necessary to account for complex underlying biological processes and the lack of ground-truth in neurosciences. We even argue that a quantitative benchmarking of methods on synthetic data without a qualitative and statistical analysis on experimental dataset make them unsuitable. Future efforts will focus on improving our simulator by modeling neuron ensemble effects, and go further in an automatic pipeline to break the neural code.

Fig. 1 :

 1 Fig. 1: Single Fluorescence trace of a neuron from Hydra Vulgaris. a) Calcium uptake of Hydra neurons displayed using Icy software [8]. b) Temporal response of one Rythmic Potential (RP) neuron with the spike pattern extracted manually.

Fig. 3 :

 3 Fig. 3: Accuracy of deconvolution methods over simulation parameters. (a) Variation of the detector gain α (i.e. level of noise). (b) Variation of the firing rate in burst. (c) Variation of the baseline amplitude.

Fig. 4 :

 4 Fig. 4: Spike Inference techniques performed on real Hydra Vulgaris' Rythmic potential (RP) and Contraction Burst (CB) neuron fluroscence traces.

Fig. 5 :

 5 Fig. 5: Despiking of the Hydra Vulgaris' dataset. a) Rasterplots of inferred spike activities show the homogeneity of the inferred spike distribution. b) Global Relaxed Jaccard Correlation.

OASIS


Marie Laviron

Maxime Petit

Eléonore Weber-Delacroix

Alexis J Combes

Arjun Rao Arkal

Sandrine Barthélémy

Tristan Courau

David A Hume

Christophe Combadière

Matthew F Krummel

Arjun Rao

Alexandre Boissonnas
email: alexandre.boissonnas@upmc.fr

Tumor-associated macrophage heterogeneity is driven by tissue territories in breast cancer

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Every adult tissue contains an abundant population of resident tissue macrophages distributed in a regular pattern with intimate interactions with other cells and extracellular matrix [START_REF] Guilliams | Establishment and Maintenance of the Macrophage Niche[END_REF][START_REF] Hume | The mononuclear phagocyte system[END_REF][START_REF] Hume | The Mononuclear Phagocyte System: The Relationship between Monocytes and Macrophages[END_REF]. Macrophages are also the most abundant immune cells of the tumor microenvironment (TME), composed of the stroma and the tumor parenchyma, and have been associated with bad prognosis in most cancers, including breast cancer [START_REF] Condeelis | Macrophages: obligate partners for tumor cell migration, invasion, and metastasis[END_REF][START_REF] Qian | Macrophage Diversity Enhances Tumor Progression and Metastasis[END_REF]. As a consequence, many proposed treatment modalities are directed towards suppression of tumor-associated macrophage (TAM) numbers or altering their function (Reviewed in [START_REF] Cassetta | Targeting macrophages: therapeutic approaches in cancer[END_REF]). Macrophage adaptation to the tumor environment has commonly been discussed in terms of a spectrum of polarization states from anti-tumor M1 to pro-tumor M2 (reviewed in [START_REF] Locati | Diversity, Mechanisms, and Significance of Macrophage Plasticity[END_REF]). However, the validity of the polarization model in cancer has been widely questioned [START_REF] Hume | Transcriptomic analysis of mononuclear phagocyte differentiation and activation[END_REF][START_REF] Mujal | Holistic Characterization of Tumor Monocyte-to-Macrophage Differentiation Integrates Distinct Immune Phenotypes in Kidney Cancer[END_REF][START_REF] Xue | Transcriptome-based network analysis reveals a spectrum model of human macrophage activation[END_REF]. Therefore, a proper understanding of TAM intrinsic diversity related to their spatial localization over time is still warranted.

During tumor growth, both tissue-resident macrophages and newly recruited monocyte-derived macrophages can give rise to TAMs which instigated interrogations of their respective role in tumor progression [START_REF] Franklin | The Cellular and Molecular Origin of Tumor-associated Macrophages[END_REF]Laviron and Boissonnas, 2019;[START_REF] Loyher | Macrophages of distinct origins contribute to tumor development in the lung[END_REF][START_REF] Müller | Resident microglia, and not peripheral macrophages, are the main source of brain tumor mononuclear cells[END_REF][START_REF] Zhu | RRID:AB_394003 Anti-mouse CD11b (M1/70) BD Biosciences Cat#563553; RRID:AB_2738276 Anti-mouse CD11b (M1/70) BD Biosciences Cat#564985; RRID:AB_2739033 Anti-mouse CD192 (SA203G11) Biolegend Cat#150605; RRID:AB_2571913 Anti-mouse CD206 (C068C2) Biolegend Cat#141723[END_REF].

Within specific tissue territories, defined as the macrophage niche [START_REF] Guilliams | Does niche competition determine the origin of tissue-resident macrophages?[END_REF], mutual repulsion and competition for available growth factors may account for the distribution of macrophages [START_REF] Guilliams | Establishment and Maintenance of the Macrophage Niche[END_REF][START_REF] Hume | The Mononuclear Phagocyte System: The Relationship between Monocytes and Macrophages[END_REF]. The link between tissue niche, ontogeny and macrophage polarization has been well described for different organs at steady state [START_REF] Bonnardel | Stellate Cells, Hepatocytes, and Endothelial Cells Imprint the Kupffer Cell Identity on Monocytes Colonizing the Liver Macrophage Niche[END_REF][START_REF] Calderon | The pancreas anatomy conditions the origin and properties of resident macrophages[END_REF]. In the mammary gland, at least two populations of resident macrophages have been described : stromal and ductal macrophages (SM and DM respectively), which are essential for the functions of the tissue [START_REF] Dawson | Tissue-resident ductal macrophages survey the mammary epithelium and facilitate tissue remodelling[END_REF][START_REF] Gouon-Evans | Requirement of macrophages and eosinophils and their cytokines/chemokines for mammary gland development[END_REF][START_REF] Jäppinen | Fetal-derived macrophages dominate in adult mammary glands[END_REF][START_REF] Wang | Tissue-resident macrophages promote extracellular matrix homeostasis in the mammary gland stroma of nulliparous mice[END_REF]. Both were shown to derive from embryonic precursors but SM were slowly replaced by monocytes over life while the majority of DM were replaced in a spike at puberty and then established as long-lived cells [START_REF] Dawson | Tissue-resident ductal macrophages survey the mammary epithelium and facilitate tissue remodelling[END_REF][START_REF] Jäppinen | Fetal-derived macrophages dominate in adult mammary glands[END_REF].

Herein, we investigated the link between spatial distribution in different territories of breast cancer and the heterogeneity of TAMs. We associated single cell RNA sequencing (scRNAseq) to spatial mapping by imaging using a multiple transgenic fluorescent reporter mouse to discriminate mammary macrophage subsets in a spontaneous and a transplanted breast tumor model [START_REF] Engelhardt | Marginating Dendritic Cells of the Tumor Microenvironment Cross-present Tumor Antigens and Stably Engage Tumor-Specific T cells[END_REF][START_REF] Guy | Induction of mammary tumors by expression of polyomavirus middle T oncogene: a transgenic mouse model for metastatic disease[END_REF]. Our analysis reveals that ductal and stromal TAMs follow distinct differentiation pathways. Further functional diversity of these TAM subsets was detected according to different anatomical structures of the stroma and tumor malignancy. DM associated with advanced tumor stage were potent phagocytes but did not show efficient capacity to activate CD8 T cells compared to SM. This spatial heterogeneity was lost in an orthotopic tumor model. Finally, similar subsets were identified in human breast tumor using a published scRNAseq dataset [START_REF] Pal | A single-cell RNA expression atlas of normal, preneoplastic and tumorigenic states in the human breast[END_REF]. Altogether, those results provide insight into the spatial heterogeneity of TAMs and the role of the local microenvironment to define their nature, which will help refining therapeutic strategies of selective TAM targeting.

Results

Breast macrophage niches evolve during tumor development

In mammary tissue, two different resident macrophage subsets were described in the literature associated with distinct locations: CD11b -ductal macrophages (DM) directly in contact with the epithelium of the mammary ducts and CD11b + stromal macrophages (SM) in the conjunctive and adipose mammary tissues [START_REF] Bijnen | Gland Macrophages: Reciprocal Control and Function within Their Niche[END_REF][START_REF] Dawson | Tissue-resident ductal macrophages survey the mammary epithelium and facilitate tissue remodelling[END_REF][START_REF] Jäppinen | Fetal-derived macrophages dominate in adult mammary glands[END_REF]. We sought to determine the evolution of those niches in the MMTV-PyMT-mCherry model (named herein PyMT) which develop spontaneous multi-focal tumors from the mammary ductal epithelium at different stages recapitulating human breast tumor development [START_REF] Attalla | Insights from transgenic mouse models of PyMTinduced breast cancer: recapitulating human breast cancer progression in vivo[END_REF][START_REF] Engelhardt | Marginating Dendritic Cells of the Tumor Microenvironment Cross-present Tumor Antigens and Stably Engage Tumor-Specific T cells[END_REF][START_REF] Lin | Progression to malignancy in the polyoma middle T oncoprotein mouse breast cancer model provides a reliable model for human diseases[END_REF]. Up to 3 different fluorescent reporters dedicated to macrophage identification were added to the PyMT background: the MacApple transgene (Csf1r mApple , [START_REF] Hawley | Csf1r -mApple Transgene Expression and Ligand Binding In Vivo Reveal Dynamics of CSF1R Expression within the Mononuclear Phagocyte System[END_REF]), the MacBlue binary transgene (ΔCsf1r-gal4/vp16/UAS-ECFP, referred to as ΔCsf1r ECFP , [START_REF] Ovchinnikov | Expression of Gal4-dependent transgenes in cells of the mononuclear phagocyte system labeled with enhanced cyan fluorescent protein using Csf1r-Gal4VP16/UAS-ECFP double-transgenic mice[END_REF]), and the knock-in Cx3cr1 EGFP/+ [START_REF] Jung | Analysis of Fractalkine Receptor CX3CR1 Function by Targeted Deletion and Green Fluorescent Protein Reporter Gene Insertion[END_REF] (Fig. 1A). These three different reporters discriminate distinct macrophage subsets according to the relative expression of Apple, ECFP and EGFP (Laviron et al., 2019). We used this fluorescent signature along with a panel of antibodies and performed tdistributed stochastic neighbor embedding (tSNE) reduction to monitor monocyte and macrophage accumulation in the mammary tissue from non-PyMT mice, 3-month-old PyMT with no detectable tumor at palpation, considered as pre-tumoral tissue and 6 month-old PyMT with palpable tumor nodules (Fig. 1B). Healthy mammary tissue (non-PyMT and pre-tumoral) contained mainly CD11b + Apple + macrophages and a small proportion of Ly6C hi ECFP + EGFP low monocytes and CD11b -EGFP + macrophages. The majority of CD11b + Apple + macrophages were also MHC-II + (IA/IE) CD206 + and CCR2 + (FigS. 1A). Two-photon imaging of fresh pre-tumoral mammary tissue confirmed that Apple + cells with typical macrophage morphology were located in the adipose and surrounding connective tissue defining a stromal population of macrophages associated to the mammary fat (SM type I) (Fig. 1C). Few round shaped ECFP + cells were detected in the connective tissue, defining a second stromal population (SM type II). CD11b -MHC-II + EGFP + macrophages corresponded to DM [START_REF] Bijnen | Gland Macrophages: Reciprocal Control and Function within Their Niche[END_REF][START_REF] Dawson | Tissue-resident ductal macrophages survey the mammary epithelium and facilitate tissue remodelling[END_REF][START_REF] Wang | Tissue-resident macrophages promote extracellular matrix homeostasis in the mammary gland stroma of nulliparous mice[END_REF]. Accordingly, these EGFP + cells were intra-epithelial (Fig. 1C, Movie 1). In tumoral tissue, CD11b -EGFP + macrophages massively accumulated along with the apparition of CD11b + EGFP + with lower expression of MHC-II (Fig. 1B, FigS. 1A). These cells did not remain only intra-epithelial but were mostly in close contact with tumor cells at all stages of tumor progression hence we still considered them as ductal TAMs (Fig. 1C, D, FigS. 1B). Based on the description of the tumor stages by Lin and colleagues, hyperplastic lesions were characterized by increased mammary epithelium density with a detectable ductal lumen [START_REF] Lin | Progression to malignancy in the polyoma middle T oncoprotein mouse breast cancer model provides a reliable model for human diseases[END_REF]. The EGFP + ductal TAMs were elongated mostly intra-epithelial, parallel to the basement membrane of the mammary duct (Fig. 1E

TAM composition shifts with tumor malignancy

EGFP + ductal TAMs comprised separable CD11b -and CD11b + subsets. To determine whether CD11b + are associated with specific location or tumor progression, we quantified by fluorescent imaging CD11b expression specifically on EGFP + cells associated to the three different stages of tumor evolution (Fig. 2A,B). CD11b expression was significantly higher in EGFP + cells associated with malignant lesions as compared to hyperplastic, consistent with the CD11b -phenotype of intra-epithelial macrophages in pre-tumoral tissue (Fig. 2B). Moreover, this increase was confirmed using flow cytometry (Fig. 2C). The relative proportion of CD11b + and CD11b -among the whole TAM compartment strikingly correlated with tumor weight but not with the age of the mouse (Fig. 2D). These results show that TAM composition shifts between hyperplastic and malignant tumor lesions.

Stromal and ductal TAMs follow distinct differentiation paths

To further characterize these different TAM subsets, three fractions of myeloid cells (total, CD11b + or CD11b -enriched) were sorted and subjected to scRNAseq. We sorted the different subsets from a pool of tumors harvested from different PyMT mice expecting to cover the different stages of tumor development (Fig. 3A). The three fractions were normalized, combined and a Louvain graphbased clustering was performed which identified 14 clusters, represented using UMAP dimensional reduction method (Fig. 3A). Differentially expressed genes were used to annotate each cluster (FigS.

2A

). First, consistent with protein expression, CD11b transcript (Itgam) was only observed in CD11b +enriched compartment (Fig. 3B,FigS. 2B,C,D). The comparison of cell distribution showed that CD11b + and CD11b --enriched fractions clustered in different compartments and that the bulk overlapped both CD11b + and CD11b --enriched fractions (FigS. 2B). Based upon a signature score from the combination of Adgre1 (F4/80), Fcgr1 (CD64) and Csf1r (CSF1R) expression, CD11b -and CD11b + compartments were both exhibiting macrophage signature (Fig. 3B,FigS. 2B,C). As expected (FigS.

1A), H2-Ab1 encoding a class II MHC gene was more highly expressed in CD11b -ductal TAMs (Fig.

3B, FigS. 2C, D).

In order to support our findings based on microscopy, we first applied to our dataset the steady state DM and SM signatures as well as lactation-and weaning-associated macrophage signatures defined by Dawson and colleagues [START_REF] Dawson | Tissue-resident ductal macrophages survey the mammary epithelium and facilitate tissue remodelling[END_REF](Fig. 3C, Table S1). DM signature mostly covered the CD11b -compartment including clusters 0, 4, 5, 8, 9, and 11, but also clusters 2 and 3 that belonged to a mix of CD11b + and CD11b -compartments. The lactation macrophage signature which reflects an active proliferation state of DM covered the same clusters than the steady state DM.

However the weaning macrophage signature reflecting DM participating in the involution of the mammary tissue through phagocytosis of milk-producing cells [START_REF] Dawson | Tissue-resident ductal macrophages survey the mammary epithelium and facilitate tissue remodelling[END_REF], covered mostly cluster 2 and a fraction of cluster 3 (Fig. 3C). SM signature mainly covered cluster 1 (Fig. 3C), but also a fraction of clusters 8 and 2 suggesting that these two clusters could be composed of both ductal and stromal macrophages or macrophages ongoing a transitory state between the stroma and the mammary duct. Dawson and colleagues also identified DC1 and DC2 signatures from sorted DC of the mammary tissue (Table S1). These signatures applied to the PyMT dataset associated cluster 7 to the DC2 signature, and cluster 10 and 13 to the DC1 signature (FigS. 2E, F). Finally, cluster 12 corresponded to a contamination of granulocytes.

Beyond the relationship established through the DM and SM signature overlaps, we next sought to match subsets identified by scRNAseq with the different TAM locations identified by imaging. We conclude that cluster 6 fits with ECFP + Ly6C high classical monocytes (Cyan circle) based on the correlation Itgam + /CD11b + , Ly6c2 + /Ly6C + , Cx3cr1 low /EGFP low (Fig. 3D) and supported by Ccr2 expression (FigS. 2F). Cluster 6 also included a small proportion of known non-classical monocyte markers Nr4a1 + , Spn + , Treml4 + (FigS. 2F). Clusters 0, 4, 5, 8, 9, 11 fitted with the EGFP + CD11b - hyperplastic ductal associated macrophages (dark green circle) based on the correspondence of expression between transcript and protein Itgam low /CD11b + , Cx3cr1 high /EGFP high , H2-Ab1 high /IA-IE high .

Clusters 2 and 3 were linked to EGFP + CD11b + malignant TAMs (light green circle) based on the correlation Itgam +/-/CD11b int , H2-Ab1 low /IA-IE low , Cx3cr1 high /EGFP high . Finally, cluster 1 was linked to the stromal type I TAMs (red circle) based on the correlation Itgam + /CD11b + , Mrc1 + /CD206, H2-Ab1 high /IA-IE high , Cx3cr1 -/EGFP -(Fig. 3D).

To address the relationship between the different macrophage subsets we performed a pseudotime trajectory analysis to model the differentiation trajectories of each subset (Fig. 3E). Franklin and colleagues used parabiosis experiment to confirm that peripheral monocytes contribute significantly to the accumulation of TAMs in the PyMT model [START_REF] Franklin | The Cellular and Molecular Origin of Tumor-associated Macrophages[END_REF]. Hence, we identified monocytes as the root node of the trajectory and defined two different paths evolving from this root. The stromal TAM path was linking monocytes to stromal TAMs and the ductal path linking monocytes to ductal TAMs through cluster 8 (Fig. 3E). Interestingly cluster 8 shared both stromal and ductal signature supporting this cluster as an intermediate state of stromal TAMs at the vicinity of the ductal epithelium.

While cluster 6 clearly identified monocytes, we were not able to clearly discriminate both subsets of stromal TAMs. To distinguish their transcriptomic profiles, we focused our analysis on cluster 1 and 6 corresponding to stromal TAMs and monocytes respectively (Fig. 3C,D) and performed a de novo clustering analysis to increase resolution. This sub-clustering generated 4 distinct clusters (FigS.

3A).

As previously, we correlated the expression of Ly6c2, Cx3cr1, H2-Ab1 and Mrc1 transcripts with their protein expression on the gated population. We recovered monocytes (Ly6c2 + /Ly6C + , H2-Ab1 low /IA-IE low ), and Apple + SM type I (Mrc1 + /CD206 + ) and we were able to assign the two other clusters to the ECFP + SM type II subset according to Cx3cr1 + /EGFP + , H2-Ab1 + /IA-IE + and Mrc1 + /CD206 + , fitting Previous studies have highlighted different markers discriminating TAM subsets in different tumor types in both human and mouse. We addressed whether these markers were also linked to different populations in our model. Hes1, Folr2 and Lyve1 expression was attributed to an embryoniclike macrophage profile [START_REF] Mulder | Cross-tissue single-cell landscape of human monocytes and macrophages in health and disease[END_REF][START_REF] Sharma | Onco-fetal Reprogramming of Endothelial Cells Drives Immunosuppressive Macrophages in Hepatocellular Carcinoma[END_REF], while Trem2+ and Spp1+ macrophages were of monocytic origin, and associated with pro-tumoral functions [START_REF] Cheng | A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells[END_REF][START_REF] Katzenelenbogen | Coupled scRNA-Seq and Intracellular Protein Activity Reveal an Immunosuppressive Role of TREM2 in[END_REF][START_REF] Mulder | Cross-tissue single-cell landscape of human monocytes and macrophages in health and disease[END_REF][START_REF] Nalio Ramos | Tissue-resident FOLR2+ macrophages associate with CD8+ T cell infiltration in human breast cancer[END_REF]. We found Folr2 and Lyve1 specifically in Mrc1+ SM; Hes1 expression was found in hyperplastic ductal TAMs, along with Cadm1.

Trem2 and Spp1 were observed in malignant TAMs (cluster 2), consistent with previous observations

(Fig. 3F, FigS. 3B).
In conclusion, TAMs accumulate through monocyte recruitment and local proliferation and progressively acquire specific transcriptomic signatures depending on their localization in distinct territories of the tumor stroma and state of tumor malignancy.

TAM niches are associated to functional diversity.

To further investigate whether the spatial heterogeneity of TAMs is associated to a specific functional profile of these different TAM subsets, we looked at the differential gene expression (DGE) between all clusters and performed a Gene Set Enrichment Analysis (GSEA)(Fig. 4A). Hyperplastic DM from cluster 4 exhibited a highly cycling score, confirming the local proliferation of CD11b -TAMs. The intermediate cluster between monocytes and DM (cluster 8) was associated with inflammatory response and strong IFN responsive scores (Fig. 4A). SM (cluster 1) were associated with a relatively active adipogenesis, fatty acid metabolism and response to hormones, which is consistent with their location around the adipocytes. Hyperplastic DM leading towards malignant DM presented specific TGFβ and TNFα signaling pathway scores and were also associated to the strongest hypoxic score. Finally, malignant DM were associated with biological hallmarks mostly linked to active metabolism (glycolysis, oxidative phosphorylation, protein secretion, fat metabolism), but also tumor-associated environment (EMT transition, apoptosis, angiogenesis, hypoxia, reactive oxygen species) (Fig. 4A). In order to better discriminate the functions among stromal TAMs, we used the DGE among the sub-clustering obtained in FigS. 3A for a second GSEA (FigS. 4A). Ly6C high monocytes represented the most active subset with high score in IFN responses, glycolysis and oxidative phosphorylation reflecting active energetic metabolism. ECFP + SM type II showed a high angiogenesis signature which might reflect their recent infiltration into the tumor.

TAM heterogeneity is commonly discussed in terms of an M1/M2 polarization spectrum. To test this model in our data, we selected a list of genes commonly related to macrophage polarization [START_REF] Biswas | Macrophage plasticity and interaction with lymphocyte subsets: cancer as a paradigm[END_REF][START_REF] Cassetta | Human Tumor-Associated Macrophage and Monocyte Transcriptional Landscapes Reveal Cancer-Specific Reprogramming, Biomarkers, and Therapeutic Targets[END_REF][START_REF] Cheng | A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells[END_REF][START_REF] Roberts | Critical Role for CD103+/CD141+ Dendritic Cells Bearing CCR7 for Tumor Antigen Trafficking and Priming of T Cell Immunity in Melanoma[END_REF]) (Table S2). To refine the analysis and seek further markers of TAM heterogeneity, we selected a set of transcripts encoding cell surface molecules involved in homing, adhesion, and migration, obtained from the HUGO database (Table S3). Hierarchical clustering of the relative expression of their transcripts confirmed that TAM subsets display specific patterns changing with tissue localization and tumor progression (Fig. 4B). We sought to identify the inherent TF involved according to the transcriptomic profile of each cluster (FigS. 4E). Putative role STAT1, STAT2, IRF7 and IRF9, associated to IFNγresponse [START_REF] Kovarik | Stat1 combines signals derived from IFNgamma and LPS receptors during macrophage activation[END_REF] were associated with the intermediate monocyte/DM cluster (cluster 8). STAT3, associated to immune suppression [START_REF] Hughes | The Multifaceted Role of STAT3 in Mammary Gland Involution and Breast Cancer[END_REF][START_REF] Hughes | Conditional deletion of Stat3 in mammary epithelium impairs the acute phase response and modulates immune cell numbers during post-lactational regression[END_REF], but also in mammary gland involution [START_REF] Sargeant | Stat3 controls cell death during mammary gland involution by regulating uptake of milk fat globules and lysosomal membrane permeabilization[END_REF], was found in the cluster of hyperplastic DM leading to malignant DM. The proliferation state of hyperplastic DM from cluster 4 was also confirmed by putative role of MYC and E2F4.

Correlation of expression between these genes across all

To gain further insights on their function beyond their different transcriptomic profile, we sorted each TAM subset from a pool of PyMT tumors loaded them with exogenous OVA peptide SIINFEKL and co-cultured them with OT-I T cells for 3 days (Fig. 4C). We measured T cell proliferation and IFNγ production as a readout of their activation. In the presence of OVA peptide SM type I were the most prone to induce T cell proliferation and IFN production by OT-I cells. (Fig. 4D). Co-culture of TAMs with sorted tumor cells did not show significant impact on tumor cell proliferation or survival (data not shown).

Taking advantage of mCherry expression by tumor cells, we next profiled the capacity of each TAM subset to participate in tumor cell clearance and thus exhibit mCherry fluorescence. Malignant CD11b + DM harbored the highest phagocytic potential among all subsets (Fig. 4E). Imaging the phagocytic EGFP + TAMs observed previously (FigS. 1C) together with CD11b staining strengthened this observation (Fig. 4F). Despite this high phagocytic ability of the mCherry OVA tumor cells, no proliferation of OT-I was detected in the absence of exogenous OVA peptide for any TAM subsets (data not shown). This is consistent with previous observations showing that highly phagocytic macrophages are not potent antigen-presenting cells [START_REF] Broz | Dissecting the Tumor Myeloid Compartment Reveals Rare Activating Antigen Presenting Cells, Critical for T cell Immunity[END_REF][START_REF] Roberts | Tissue-Resident Macrophages Are Locally Programmed for Silent Clearance of Apoptotic Cells[END_REF]. Overall, we conclude that the localization of TAM subsets drives their specific transcriptional programs and functional capacities in the tumor microenvironment.

TAM heterogeneity is reduced in orthotopic model

We next addressed whether similar TAM heterogeneity can be recovered in an orthotopic tumor model. We performed orthotopic injection of breast tumor cells into the mammary fat pad as a model of the malignant stage of tumor development. PyMT cell line derived from a spontaneous MMTV-PyMT-mCherry tumor were injected to MacApple x MacBlue x Cx3cr1 EGFP/+ mice and tumor growth and TAM composition were monitored (Fig. 5A). Tumor nodules were palpable after 11 days and then grew exponentially to reach 4000 mm 3 within a further 10 days (Fig. 5B). Among monocytes only ECFP + EGFP low signature was recovered whereas among macrophages the three fluorescent signatures, also observed in the spontaneous tumors, were present (Fig. 5C). No CD11b -EGFP + cells corresponding to DM were recovered in this model and ECFP + TAMs and EGFP + TAMs represented the most abundant subsets with a progressive accumulation with tumor growth (Fig. 5D). The Apple + TAMs remained the less abundant subset and did not significantly accumulate over tumor growth (Fig. 5D).

CD206 did not discriminate a specific subset of stromal TAMs, as observed in the spontaneous model, as the expression was similar in all subsets (FigS. 5A). We analyzed the spatial distribution of these different TAM subsets by histological analysis. The structure of the tumor appeared as a dense aggregate of tumor cells with no distinction of stroma or hyperplastic mammary epithelium as observed in the spontaneous PyMT model (Fig. 5E). The tumor was homogeneously infiltrated by the three fluorescent cell subsets.

We again performed scRNAseq on sorted myeloid cells (Fig. 5F). Louvain graph-based clustering identified 5 clusters (Fig. 5G). Clusters 3, 4, 5 exhibited DC-associated signatures as defined by Dawson et al. (Fig. 5G,H, FigS. 5B) [START_REF] Dawson | Tissue-resident ductal macrophages survey the mammary epithelium and facilitate tissue remodelling[END_REF]. Cluster 1 corresponded to Ly6C hi monocytes (Itgam + , Ly6c2 + , Ccr2 + ), while cluster 2 showed an intermediate monocyte/macrophage signature. Cluster 0 was the only cluster representing TAMs. DM and SM signatures both covered this cluster (Fig. 5H) suggesting that only one main TAM differentiation path occurs in the orthotopic model in accordance with the observation that all TAMs develop in a rather homogeneous TME, mainly composed of tumor cells in this model. Additionally, the expression of Trem2, Cadm1, Folr2 and Mrc1 was recovered in this TAM cluster, although their level of expression seemed to identify different cells within the cluster (Fig. 5I). The GSEA based on the DGE of the 5 clusters revealed that cluster 0 was highly inflammatory and metabolically active, with a proliferative signature, resembling the cluster associated to malignant TAMs in the spontaneous model (cluster 2) (FigS. 5C). Consistently, the projection of the signatures extracted from clusters 0, 1 and 2 of the orthotopic model to the UMAP of the spontaneous PyMT model confirmed the similarities with the cluster of malignant TAMs and hyperplastic TAMs (FigS. 5D). We conclude that orthotopic models negate the impact of specific tissue location and drive TAMs toward a functional signature representative of advanced tumor malignancy.

Human breast TAMs exhibit similar diversity

As the MMTV-PyMT model has been described to recapitulate the different stages of development of human breast cancer, we sought to determine whether the TAM heterogeneity identified in this model could be reflected in human tumors. We used scRNAseq data from Visvader's group containing a diversity of cell types sequenced from different breast tumors [START_REF] Pal | A single-cell RNA expression atlas of normal, preneoplastic and tumorigenic states in the human breast[END_REF]. We focused our analysis on monocytes and macrophages by selecting them based on CSF1R, CD68 or FCGR1A expression and generated a UMAP representation of the concatenated tumor-associated myeloid cells.

Louvain graph-based clustering identified 10 clusters (Fig. 6A). We applied the combination of markers that were used to identify distinct TAM subsets in PyMT mouse to the human dataset (Fig. 6A). While those signatures were in general more heterogenous, we retrieved a similar association. Monocytes were identified by expression of ITGAM, CD14 and CCR2. FOLR2 + cells co-expressed LYVE1 and MRC1 as observed in the stromal/adipose TAMs in the PyMT model (Fig. 6A,B). A population of ITGAM low HES1 + CADM1 + was recovered, suggesting that a profile similar to Hyperplastic DM in PyMT mice can be found in human breast tumors. Finally, advanced TAMs were identified by TREM2, SPP1 and MARCO expression. An intermediate population expressing markers from both hyperplastic-like and advanced-like TAMs, was identified (Fig. 6A,B). Similar observations were recently made by Nalio Ramos and colleagues showing that FOLR2 identify macrophages in the stroma in both healthy mammary tissue and malignant breast cancer while TREM2 and SPP1 expression is specific to TAMs localized in the tumor nest [START_REF] Nalio Ramos | Tissue-resident FOLR2+ macrophages associate with CD8+ T cell infiltration in human breast cancer[END_REF]. Stromal macrophages exhibited an androgen response-associated transcriptomic signature which is coherent with their association to the adipose stroma of the mammary tissue (Fig. 6C). As in mice, an IFN-responsive cluster was also recovered in human tumors, corresponding to advanced-like TAMs. Altogether, these observations suggest that similar niche-associated patterns are found in human breast tumors.

Discussion

Spontaneous tumors from the PyMT model generate multiple tumors from the mammary epithelium representative of the different stages of breast cancer in human tumors [START_REF] Attalla | Insights from transgenic mouse models of PyMTinduced breast cancer: recapitulating human breast cancer progression in vivo[END_REF][START_REF] Lin | Progression to malignancy in the polyoma middle T oncoprotein mouse breast cancer model provides a reliable model for human diseases[END_REF]. Hence, this model provides both the spatial diversity of the mammary tissue and the temporal evolution of tumor transformation. Using three different reporter transgenes, we were able to distinguish multiple subsets of TAMs, associated to distinct spatial and temporal territories. These different TAMs harbored conserved transcriptomic signatures of the two main macrophage populations isolated from the healthy mammary tissue showing that environmental perturbation induced by tumor development does not fully reprogram TAMs from their original function and that TAM heterogeneity is spatially poised.

CD11b -ductal TAMs, mainly detected in early tumors were linked to monocytes through a differentiation path associated to IFN signaling signatures, known to be involved in macrophage differentiation and response to immune checkpoint blockade therapy [START_REF] Benci | Tumor Interferon Signaling Regulates a Multigenic Resistance Program to Immune Checkpoint Blockade[END_REF][START_REF] Delneste | Interferon-gamma switches monocyte differentiation from dendritic cells to macrophages[END_REF]. This subset also exhibited lactation macrophage signature, consistent with the idea that hyperplastic lesion resembles mammary tissue in active proliferation in which DM and epithelial cells exert mutual benefits [START_REF] Dawson | Tissue-resident ductal macrophages survey the mammary epithelium and facilitate tissue remodelling[END_REF][START_REF] Van Nguyen | Colony stimulating factor-1 is required to recruit macrophages into the mammary gland to facilitate mammary ductal outgrowth[END_REF]. Accordingly, the active proliferation of this subset suggests that they partially accumulate through local proliferation. We identified a secondary path leading to malignant ductal TAMs with high Trem2 expression. Accordingly, TREM2 has been associated with bad prognosis and was suggested to regulate the immunosuppressive functions of TAMs of monocytic origin [START_REF] Katzenelenbogen | Coupled scRNA-Seq and Intracellular Protein Activity Reveal an Immunosuppressive Role of TREM2 in[END_REF][START_REF] Mulder | Cross-tissue single-cell landscape of human monocytes and macrophages in health and disease[END_REF]. Although this subset reflects the presence of more advanced tumor lesions, these TAMs do not necessarily exert only pro-tumoral function. For instance, they harbored the most potent phagocytic capacity and similar signature to DM involved in tissue involution after weaning. TNF-and IFN-responsive clusters were identified in both human and mouse tumors suggesting that these molecular axes are important paths in TAM activation and could represent promising therapeutic targets. Additionally, both pathways were recently shown to a define a specific tumor immune archetype [START_REF] Combes | Discovering dominant tumor immune archetypes in a pan-cancer census[END_REF].

A third independent path linking monocytes to stromal TAMs was identified. Stromal type I TAMs harbored potent T cell activation capacity in vitro, although they phagocyted less tumor cells than the CD11b + ductal TAMs. This is likely to be associated with their respective tissue localization with stromal type I TAMs restricted to the border of the tumor and in adipose islets while CD11b + ductal TAMs are found in advanced tumor lesions. Similar observations were recently made in human breast cancer where FOLR2 + TAMs were positioned near blood vessels at the tumor periphery and co-localized with CD8 T cells [START_REF] Nalio Ramos | Tissue-resident FOLR2+ macrophages associate with CD8+ T cell infiltration in human breast cancer[END_REF]. We confirmed that the stromal type I TAMs exhibited higher capacity to activate OT-I T cells in vitro compared to other subsets only when adding exogenous cognate antigen. However OT-I adoptive transfer in vivo has been shown to ineffective at rejecting the tumor [START_REF] Engelhardt | Marginating Dendritic Cells of the Tumor Microenvironment Cross-present Tumor Antigens and Stably Engage Tumor-Specific T cells[END_REF].

scRNA-seq has some limitations as a method to detect genuine heterogeneity in populations.

Even the most highly-expressed transcripts are detected in only a subset of cells [START_REF] Summers | Network analysis of transcriptomic diversity amongst resident tissue macrophages and dendritic cells in the mouse mononuclear phagocyte system[END_REF] and clustering of cell types depends heavily on assumptions about the data structure. However, our clustering was mostly confirmed at the protein level by flow cytometry and recovered the main subsets observed by imaging. We showed that most TAM transcriptomic signatures poorly recovered the M1/M2 dichotomy, as shown recently in other mouse models [START_REF] Mujal | Holistic Characterization of Tumor Monocyte-to-Macrophage Differentiation Integrates Distinct Immune Phenotypes in Kidney Cancer[END_REF]. CD206 is usually associated to M2-like phenotype but is widely expressed by subsets of resident tissue macrophages including those of adipose tissue [START_REF] Arendt | Obesity promotes breast cancer by CCL2-mediated macrophage recruitment and angiogenesis[END_REF][START_REF] Silva | Vasculature-associated fat macrophages readily adapt to inflammatory and metabolic challenges[END_REF][START_REF] Summers | Network analysis of transcriptomic diversity amongst resident tissue macrophages and dendritic cells in the mouse mononuclear phagocyte system[END_REF][START_REF] Wentworth | Pro-inflammatory CD11c+CD206+ adipose tissue macrophages are associated with insulin resistance in human obesity[END_REF].

Although those TAMs were not in close contact with tumor nodules, adipose-associated macrophages are known to produce several inflammatory cytokines promoting tumor development [START_REF] Arendt | Obesity promotes breast cancer by CCL2-mediated macrophage recruitment and angiogenesis[END_REF][START_REF] Faria | Obesity and Breast Cancer: The Role of Crown-Like Structures in Breast Adipose Tissue in Tumor Progression, Prognosis, and Therapy[END_REF][START_REF] Picon-Ruiz | Obesity and adverse breast cancer risk and outcome: Mechanistic insights and strategies for intervention[END_REF]. In the orthotopic model CD206 is expressed in a high frequency of TAMs regardless of their localization but in human tumors, Mrc1 + stromal TAMs were present in a proportion similar to the spontaneous model. This shows that orthotopic models tend to artificially bias TAM phenotype. Indeed, TAM diversity was more restricted in the orthotopic model which we argue to be linked to the lack of niche variety and to a model representative of very progressive tumor. In the human dataset, we were able to recover a similar TAM diversity including the HES1 + , FOLR2 + and TREM2 + main subsets with functional hallmarks comparable to the ones present in the mouse PyMT tumor. These subsets were recently identified in human breast tumors to reside in distinct niches [START_REF] Nalio Ramos | Tissue-resident FOLR2+ macrophages associate with CD8+ T cell infiltration in human breast cancer[END_REF]. Those observations emphasize the importance to consider TAMs in human tumors according to their tissue localization [START_REF] Wu | A single-cell and spatially resolved atlas of human breast cancers[END_REF]. The development of mass cytometry imaging and spatial sequencing will further address this aspect.

Taken together, our data suggest that heterogeneity derived not from any specific activation state, but from the temporal profile adaptation of monocyte and macrophage to specific locations with the growing and developing tumor. Overall, our study propose that TAM heterogeneity is directly related to the spatial diversity which is shaped by the original structure of the tissue, the histological type of the tumor and its developmental grade. Therefore, it will be key in the future to consider TAM spatial distribution to develop more tailored therapeutic strategies to precisely shape the tumor immune response.

Limitations of study

Although we were able to link spatial heterogeneity of TAMs to their transcriptomic signature, the fluorescent profiles did not fully recover the granularity of TAMs reached by scRNAseq but at least was associated to the main tissue territories. The reason why TAM subsets differentially expressed the fluorescent reporters as a result of their localization is still unclear. Further functional studies will help to determine how the surrounding signals regulate fluorescent reporter expression. Animales, Environnement with the number #16890. Sample sizes were chosen to assure reproducibility of the experiments and according to the 3 Rs of animal ethic regulation. For scRNAseq, all mice were handled in accordance with NIH and American Association of Laboratory Animal Care standards, and experiments were approved by the Institutional Animal Care and Use Committee of the University of California, San Francisco.

PyMT Tumor Cell line

The PyMT cell line was derived from primary PyMT-ChOVA female breast tumor as described (You et al.). Cells were cultured in RPMI-1640 with 10% fetal bovine serum, Na-Pyruvate, antibiotic, and antimycotic (GIBCO) and incubated at 37°C in 5% CO2. For orthotopic breast cancer model, 2.5 × 10 5 cells were injected in the mammary fat pad of 20-weeks old female.

METHOD DETAILS

Mouse Tissue Digestion and Flow Staining

Healthy mammary tissue, orthotopic and spontaneous tumors were processed the same way. Mammary tissues were harvested and enzymatically digested with 0.1 mg/ml DNase I (Sigma-Aldrich), 1.5 U/ml Dispase II (Gibco), and 200 U/ml Collagenase Type IV (Gibco) in RPMI+glutamax (Gibco) for 30 minutes at 37°C under agitation. Samples were filtered using 70µm cell strainer and washed with FACS buffer (PBS, 0.5% FCS, 0.01% azide, 2mM EDTA). Cells were then washed and non-specific binding was blocked with Fcblock (BD biosciences). Cell surface proteins were then stained at 4°C for 30 minutes in Brilliant stain buffer (BD biosciences). Anti-CD45 (30-F11), anti-CD11b (M1/70), anti-Ly6C (AL21), anti-Ly6G (1A8), anti-SiglecF (E50-2440), anti-I-A/I-E (M5/114.15.2), anti-F4/80 (T45-2342) were purchased from BD Biosciences. Anti-CD64 (X54-5/7.1), anti-CD192 (SA203G11) and anti-CD206 (C068C2) were purchased from Biolegend. Cells were then washed once in FACS buffer and analyzed directly by flow cytometry. For cell cycle analysis, cells were fixed in 1% PFA for 10min at 4°C then washed in PBS, 2%FCS, 0.01% azide, 0.1% saponin. Cells were stained 1h at RT in 500μL PBS containing 2 μM topro-3 iodide (Molecular probes) and 10 μM RNase A (Boehringer). Calculation of absolute cell number was performed by adding to each vial a fixed number (10.000) of non-fluorescent 10μm polybead carboxylate microspheres (Polysciences, Niles, IL, USA) according to the formula: Nb of cells = (Nb of acquired cells x 10.000)/(Nb of acquired beads). Number of cells obtained for each sample was normalized per mg of tissue. Flow cytometry acquisition was performed on the flow cytometer FACS LSRFortessa X-20 (BD Biosciences, Franklin Lakes, NJ, USA) with DIVA (BD Biosciences) Flow Cytometry software. Flow cytometry data were analyzed with FlowJo software (Tree Star Inc., Ashland, OR, USA) or OMIQ (Omiq inc.) data analysis software for opt-SNE representations (https ://omiq.ai, Santa clara CA).

Single cell RNA sequencing

Myeloid cells were sorted using a BD FACSAria Fusion (BD Biosciences, Franklin Lakes, NJ, USA), according to alive CD45 + CD90.2 -B220 -NK1.1 -Ly6G -CD11b + or CD11c + from a mix of different tumors isolated from a pool of 4 PyMT-ChOVA mice. Among these, CD11b + and CD11b -myeloid compartment were individually sorted and processed for independent sequencing along with the bulk compartment.

For the orthotopic model, only cells from the bulk gating were sorted from a pool of 5 different mice 15 days after tumor inoculation. After sorting, all cell preparations were pelleted and resuspended at 1x10 3 cells/ml in 0.04%BSA/PBS and loaded onto the Chromium Controller (10X Genomics). Samples were processed for single-cell encapsulation and cDNA library generation using the Chromium Single Cell 3'v2 Reagent Kits (10X Genomics) according to manufacturer's protocol. The library was subsequently sequenced on an Illumina HiSeq4000 (Illumina).

Single Cell Data Processing

Sequencing data was processed using 10X Genomics Cell Ranger V1.2 pipeline. The Cell Ranger subroutine mkfastq converted raw, Illumina bcl files to fastqs. Reads were then aligned with count, using the STAR aligner [START_REF] Dobin | STAR: ultrafast universal RNA-seq aligner[END_REF] against the mm10 mouse genome. After filtering reads with redundant unique molecular identifiers (UMI), count generated a final gene-cellular barcode matrix. Both mkfastq and count were run with default parameters.

Cellular Identification and Clustering

For each sample, the gene/barcode matrix was passed to the R software package Seurat (v3.1.5) [START_REF] Stuart | Comprehensive Integration of Single-Cell Data[END_REF] for all downstream analyses. We then filtered on cells that expressed a minimum of 200 genes and required that all genes be expressed in at least 3 cells. We also removed cells that contained > 9% reads associated with cell cycle genes [START_REF] Kowalczyk | Single-cell RNA-seq reveals changes in cell cycle and differentiation programs upon aging of hematopoietic stem cells[END_REF][START_REF] Macosko | Highly Parallel Genome-wide Expression Profiling of Individual Cells Using Nanoliter Droplets[END_REF].

Count data was then log2 transformed and scaled using each cell's proportion of cell cycle genes as a nuisance factor (implemented in Seurat ScaleData function) to correct for any remaining cell cycle effect in downstream clustering and differential expression analyses. For the spontaneous PyMT model, Bulk, CD11b + and CD11b -enriched myeloid compartment data sets were merged for downstream analysis.

Principal component (PC) analysis was performed on a set of highly variable genes defined by Seurat FindVariableGenes function. Genes associated with the resulting PCs were then used for graph-based cluster identification and subsequent dimensionality reduction using UMAP. Cluster-based marker identification and differential gene expression (DGE) were performed using Seurat FindAllMarkers for all between-cluster comparisons. Few clusters enriched for Epcam-Krt identified as tumor cells and enriched in ribosomal transcripts considered as dying cells were discarded from the analysis. Subclustering was performed by selecting cells from clusters of interest and the whole same processing was applied.

Gene signature overlay

Stromal macrophage, ductal macrophage, weaning, lactation, DC1 and DC2 were downloaded from [START_REF] Dawson | Tissue-resident ductal macrophages survey the mammary epithelium and facilitate tissue remodelling[END_REF] (Table S1). We created a dataframe of gene values per cell from the input Seurat object. The calculate score was visualized on the UMAP with a saturation score between 0.5 and 0.7 for each signature. Density plot of genes and signatures were generated using the Nebulosa package (Alquicira-Hernandez and Powell, 2021).

Pseudotime analysis

Trajectory analysis was performed using the Monocle3 package (v1.0.0) [START_REF] Cao | The single-cell transcriptional landscape of mammalian organogenesis[END_REF][START_REF] Qiu | Reversed graph embedding resolves complex single-cell trajectories[END_REF][START_REF] Trapnell | The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells[END_REF]. Resolution was determined according to the clustering obtained with the Seurat package. Root node was assigned to monocytes for each analysis.

GSEA Hallmark pathway analysis

The signature of each cluster was computed into GSEA (https://www.gsea-msigdb.org) for top100 hallmark pathways with FDR q-value < 0.05. From the generated list, only the hallmark with a k/K value > 0.04 were selected. A score of 0 was set when the pathway was not found in the cluster. Heatmap was generated using the pheatmap package. k/K values of each pathway were used for the heatmap.

Hierarchical clustering of centered and scaled gene expression was performed using complete linkage and Euclidean distance.

M1/M2 gene correlation

The list of commonly M1 or M2-associated genes was obtained from [START_REF] Biswas | Macrophage plasticity and interaction with lymphocyte subsets: cancer as a paradigm[END_REF][START_REF] Cassetta | Human Tumor-Associated Macrophage and Monocyte Transcriptional Landscapes Reveal Cancer-Specific Reprogramming, Biomarkers, and Therapeutic Targets[END_REF][START_REF] Cheng | A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells[END_REF][START_REF] Roberts | Critical Role for CD103+/CD141+ Dendritic Cells Bearing CCR7 for Tumor Antigen Trafficking and Priming of T Cell Immunity in Melanoma[END_REF] (Table S2). From this complete gene list (69 genes), 61 were recovered from the dataset and 50 out of the 61 were conserved using the Seurat FindVariableGenes function. Genes with lowest sct.mean (<0.02) for the dataset were discarded. We used the R package Stats (version 4.0.3) to calculate the correlation scores for each gene in individual cells. The cor function was used to compute gene correlation score with "spearman" as method.

Resulting heatmap was generated with pheatmap. The same list was used to identify specific pattern in each cluster. The heatmap was generated using Doheatmap.

Homing molecule heatmap

Homing molecules were extracted from the HUGO database based on the manual selected lists of Regulators of G protein signaling, integrins, selectins, and chemokine receptor families (Table S3). The most differentially expressed genes among clusters from this list were determined using FindAllMarkers from Seurat and the heatmap was generated using Doheatmap.

Transcription factor analysis

We applied murine version of DoRothEA (version 1.2.1) package in combination with VIPER package (version 1.24.0) on single sample matrices (genes in rows and single cells in columns) containing normalized gene expression scores scaled gene-wise. Transcription factors with lowest confidence levels listed in the corresponding gene expression matrix were discarded from the analysis. VIPER provides a normalized enrichment score (NES) for each TF which we consider as a metric for the activity.

Resulting heatmap was generated with pheatmap.

Human data processing

scRNAseq dataset were obtained from [START_REF] Pal | A single-cell RNA expression atlas of normal, preneoplastic and tumorigenic states in the human breast[END_REF])(GSE161529). All 32 files corresponding to female tumor samples of all types were pre-processed with the same pipeline than mouse scRNAseq.

After cleaning the data, myeloid cells identified by the expression of PTPRC, CD68, CSF1R, FCGR1A of each individual tumor were subsampled. Tumors containing less than 300 myeloid cells were excluded for downstream analyses. Samples were combined using anchor-based integration method from Seurat. The integrated data were used for downstream analysis, for which the same pipeline than mouse analysis was applied.

Co-culture experiment

Each TAM subset was sorted using a BD FACS Aria III according to alive CD45 + Ly6G -Ly6C -CD64 + , CD11b level and fluorescent reporter expression from pooled tumors of PyMTFluo-ChOVA mice. OT-I T cells were isolated from axillary and brachial lymph nodes and incubated in 1mM CFSE for 5min at 37°C and washed in PBS/FCS. Cells were plated at a 5:1 T cell/TAM ratio in 96 round well plate in RPMI with 10% FCS with antibiotics/antimycotics. 1µM ovalbumin peptide SIINFEKL was added to each well.

Plates were incubated for 72 hours at 37°C -5%C CO2. For the last three hours of culture, 5µg/ml of brefeldin A was added to each well for IFNγ staining.

Multi-photon imaging

The two-photon laser-scanning microscopy (TPL SM) set-up used was a 7MP (Carl Zeiss) coupled to a Ti: Sapphire Crystal multiphoton laser (ChameleonU, Coherent), which provides 140-fs pulses of nearinfrared light, selectively tunable between 680 and 1050 nm and an optical parametric oscillator (OPO-MPX, Coherent) selectively tunable between 1,050 and 1,600 nm. The NLO and the OPO beams were spatially aligned and temporally synchronized using a delay line (Coherent) allowing CARS imaging approach [START_REF] Boissonnas | Imaging resident and recruited macrophage contribution to Wallerian degeneration[END_REF]. The excitation wavelength was 840 nm for the NLO beam and 1104 nm for the OPO beam to detect the vibrational signature of lipid rich structures at a frequency of 2846 cm -1 with an emission wavelength at 678 nm. The system included a set of external nondescanned detectors in reflection with a combination of a LP-600-nm followed by LP-462-nm and LP-500-nm dichroic mirrors to split the light and collect the second harmonic generation signal (SHG) with a 417-/60-nm emission filter, ECFP with a 480-/40-nm emission filter, EGFP with a 525-/50-nm emission filter, mCherry or Apple with a 624-/40-nm emission filter and CARS or blue evans signal with a LP 645nm emission filter. For live imaging, mice were anaesthetized and maintained during the imaging period with 2% isoflurane in medical air. An incision of the skin was performed at the level of the breast tumor and an imaging window was positioned to stabilize the tissue from breathing and drifting artefacts. Local temperature was monitored and maintained at 32°C using an incubation chamber. To define the tumor vasculature, Evans Blue was injected i.v. before the imaging session. Real time movies were performed by imaging every 30s by 5 consecutive 3μm z spacing image stack (total 12μm thickness). Static 3D images were performed on 100µm thick cryosections by 1μm z spacing image stack (up to 20μm thickness). For all images the objective was a water immersion, plan apochromat ×20 (numerical aperture = 1). 3D reconstructions and mask rendering were done using Imaris software (Bitplane).

Epifluorescent imaging of tumor sections

In brief, mammary tissue and tumors were harvested and fixed in 2% PFA for 6 h and then incubated in 30% sucrose-PBS overnight at 4°C before being embedded in OCT-freezing medium (Microm microtech) and frozen at -80°C. Sectioning was completed on a HM550 Cryostat (Thermo Fisher) at -20°C; 5-μm or 100-μm sections were collected on Superfrost Plus Slides (Thermo Fisher Scientific) and stored at -20°C until use. Tissue sections were rehydrated with 0.5% Triton-PBS during 10 min.

For CD11b staining, a first block step was performed with 3% BSA solution during 2h, followed by 2h incubation at room temperature with the primary antibodies (Rat anti-mouse CD11b, clone M1/70; BD Biosciences). Secondary antibody incubation was performed using anti-rat AF647-conjugated antibody 1h30 at room temperature (Jackson ImmunoResearch Laboratories). Slides were counterstained and mounted with Vectashield Mounting Medium with DAPI (Vector Laboratories). Epifluoresent images were acquired on a Zeiss Axio Z1 fluorescent microscope (Carl Zeiss, Germany) using Zen software. adenoma, malignant, according to their histological characteristics [START_REF] Lin | Progression to malignancy in the polyoma middle T oncoprotein mouse breast cancer model provides a reliable model for human diseases[END_REF]). Acquisition and analysis settings were identical for both isotype and CD11b staining. Between three and six fields were chosen for each lesion type per mouse (n=6 mice). All histological quantifications are presented as a mean of the different fields for each mouse.

  , FigS. 1C). Adenoma lesions were characterized by tumor nodules with no more ductal lumen but a clear distinction of epithelial basement in which EGFP + ductal TAMs displayed stellate morphology at the periphery. Malignant lesions were characterized by a disrupted epithelial basement with invasive tumor organization where EGFP + ductal TAMs surrounded the tumor or were completely trapped within the tumor parenchyma (Fig. 1E, FigS. 1C). Highly phagocytic EGFP + ductal TAMs were observed located in the lumen of the mammary duct similar to the ones observed during mammary epithelium involution post weaning (FigS. 1C)(Dawson et al., 2020). Upon tumor development another macrophage subset defined as CD11b + ECFP + EGFP low (for simplicity called ECFP + ) (Fig. 1B), appeared and was mostly detected in stromal regions and less in close contact with tumor cells (Fig. 1C, E, FigS. 1B, D). These ECFP + cells were composed of both monocytes and macrophages (Fig. 1B, C). Numerous ECFP + cells were detected rolling and patrolling in the vasculature as expected for monocytes (Movie 2). They also accumulated nearby vascular regions within the tumor parenchyma and displayed a higher motility pattern compared to EGFP + cells (Movie 3) but exhibited decreased velocity as they were moving closer to tumor cells (Movie 4, FigS. 1E). We assumed that ECFP + cells represent the most recently infiltrating monocytes and macrophages that progressively migrate and differentiate in the different regions of the tumor. Finally, upon tumor expansion, the Apple + SM type I were still mostly located to the adipose tissue at the periphery of the tumor or in smaller adipose islets (Fig. 1C, FigS. 1B). To summarize, EGFP + cells were considered as ductal TAMs further distinguished by the emergence of CD11b + subset. Apple + cells represented a SM type I population associated to the mammary fat. ECFP + cells represented both monocytes and a distinct tumor-induced SM type II population. In conclusion, we were able to identify in mammary tumors several subsets of TAMs exhibiting distinct localization related to the different macrophage niches existing prior to tumor development.

  with the phenotype observed by flow cytometry (FigS. 3A, B, C). Pseudotime analysis on the subclustering of the stromal TAMs suggested a linear path from monocytes to stromal type I and stromal type II TAMs suggesting that monocytes may differentiate into both type of stromal TAMs (FigS. 3D). Following the ductal path, a circular relationship was shown in the region of hyperplastic ductal TAMs suggesting a contribution of local proliferation. This was supported by cell cycle analysis by flow cytometry using Topro3 to stain the DNA, showing that cycling CD11b -TAMs represented the vast majority of cycling macrophages (Fig. 3E, FigS. 3E). In contrast, local proliferation weakly contributed to the maintenance of SM type I macrophages which might explain their relative loss along tumor evolution (FigS. 3E). Finally, from the hyperplastic ductal TAM clusters, an alternative path throughout cluster 3 led to malignant tumor TAMs identified in cluster 2 (Fig. 3E).

  individual cells was determined. Overall, correlation scores were low and hierarchical clustering did not discriminate convincing M1-or M2associated sets of genes, with the exception of two modules (correlation score >0.7) (FigS. 4B). The first module defined by Lyve1, Il10, Cd163 and Ccl24 identified a fraction of the Mrc1 + cluster 1 corresponding to SM (FigS. 4C). The second module defined by Arg1, Arg2 and Mmp9 highlighted a fraction of malignant TAMs (FigS. 4C). Individual gene expression across the different clusters confirmed the absence of a clear-cut canonical M1 or M2 profile (FigS. 4D) suggesting that the spatial localization of TAMs does not seem to reconcile the M1 or M2 nomenclature.

Figure 1 .

 1 Figure 1. Breast macrophage niches evolve during tumor development. A. MMTV-PyMT-mChOVA were crossed with Cs1fr-mApple, ΔCsf1r-ECFP and Cx3cr1 kin/GFP to generate a multi-fluorescent macrophage reporter mouse (PyMTFluo). B. Representative opt-SNE of flow cytometry analysis of the mammary monocyte/macrophage compartment shows the expression of the different markers and fluorescent reporters over tumor development. Non-PyMT represents non tumoral mammary tissue from female fluorescent mice. Pre-tumoral tissue represents mammary tissue from < 3-month-old tumor-free female PyMTFluo mouse. Tumoral tissue represents mammary tissue from > 6-month-old tumor-bearing female PyMTFluo mice. The different macrophage subsets are colorcoded based on their relative expression of ECFP, EGFP, Apple, CD11b and Ly6C. Histogram represents the percentages of each subset of fluorescent macrophages (n= 5-8 mice per group, mean ± SD are indicated). Two-way ANOVA with Turkey post test was performed. C. TPLSM images of the mammary tissue show the distribution of the different macrophage subsets in pre-tumoral and tumoral mammary tissue from PyMTFluo mice. Collagen was detected by second harmonic generation (SHG). Adipocytes were detected by coherent anti-stokes raman spectroscopy (CARS) (Boissonnas et al., 2020). Scale bars 20μ. D. TPLSM image of mammary tumor border shows the restriction of stromal Apple + cells to adipose stroma. Scale bar 20μ. E. 3D TPLSM reconstructions show accumulation of EGFP + and ECFP + TAMs in the different types of tumor lesions; hyperplasia, adenoma, malignant ductal, and stroma. Scale bar 20μ. All panels all representative of at least n = 4 mice. See also Figure S1.

Figure 2 .

 2 Figure 2. Ductal TAMs express CD11b + in malignant tumor lesions.

Figure 3 .

 3 Figure 3. Stromal and ductal TAMs harbor localization-dependent transcriptomic profiles and differentiation paths

Figure 4 .

 4 Figure 4. TAM subsets exhibit distinct polarization

Figure 5 .

 5 Figure 5. TAM heterogeneity is associated to niche diversity

Figure 6 .

 6 Figure 6. Human breast TAMs exhibit a niche-associated TAM diversity

  ECFP, EGFP, DAPI, mCherry, and AF647 signals were acquired using a combination of LED Colibri 7 and emission filters: LED 385 nm, EmBP 450/40 for DAPI; LED 430nm, EmBP 480/40 for ECFP; LED 475 nm, EmBP 525/50 for EGFP; LED 511nm, EmBP 605/70 for Apple; LED 555nm, EmBP 647/57 for mCherry; LED 630nm, EmBP 690/50 for AF674. Image quantifications were performed using histocytometry based approach. Briefly, after masking the cells as individual objects based on DAPI, EGFP or mCherry signal for tumor nodule contouring using Imaris software, fluorescent signals from all channels and spatial coordinates were extracted and transformed into fcs file for cytometry-based analysis using FlowJo. Distance of EGFP + and ECFP + TAMs to the tumor was calculated using the distance transformation module from Imaris software. CD11b quantification was performed by surfacing TAMs based on the EGFP channel. CD11b expression on EGFP + TAMs was measured using the sum intensity of the CD11b channel parameter gated on different regions of the tumor (hyperplastic,

Acknowledgments

This work was supported by Association pour la Recherche contre le Cancer ARC "Projet labellisé", by Cancéropole Ile-de-France "programme Emergence 2020". INCa (2018-1-PLBIO-06-1) supported research activity and M.P. funding. M.L. was funded by Sorbonne University program "Idex SUPER" and Fondation pour la Recherche Médicale FRM. The authors thank UMS28 for the animal core facility and the two-photon microscopy, and Doriane Foret and Thomas Van Hecke for animal breeding. Research in the D.A.H laboratory was supported by the Mater Foundation. AJ. C, TC and MF. K were supported by UCSF Bakar ImmunoX and NIH (R01 U01).

Conflict of interest

The authors declare no competing interests.

STAR methods

RESOURCE AVAILABILITY

Lead contact

Further information and requests for resources and reagents should be directed to Alexandre Boissonnas (alexandre.boissonnas@upmc.fr).

Materials availability

This study did not generate new unique reagents.

Data and code availability

 Single-cell RNAseq data have been deposited at Gene Expression Omnibus and are publicly available as of the date of publication at the following accession number: GSE184096. GEO accession number is also listed in the key resources table.

 This paper does not report original code.

 Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

EXPERIMENTAL MODELS AND SUBJECTS DETAILS

Mice

Female MMTV-PyMT-P2A-mCherry-P2A-OVA (PyMT-ChOVA) [START_REF] Engelhardt | Marginating Dendritic Cells of the Tumor Microenvironment Cross-present Tumor Antigens and Stably Engage Tumor-Specific T cells[END_REF] for the scRNAseq were bred under specific pathogen-free conditions at the University of California, San Francisco Animal Barrier Facility. For flow cytometry and imaging analyses, Cx3cr1 EGFP/+ [START_REF] Jung | Analysis of Fractalkine Receptor CX3CR1 Function by Targeted Deletion and Green Fluorescent Protein Reporter Gene Insertion[END_REF], Csf1r-Gal4VP16/UAS-ECFP (MacBlue; [START_REF] Ovchinnikov | Expression of Gal4-dependent transgenes in cells of the mononuclear phagocyte system labeled with enhanced cyan fluorescent protein using Csf1r-Gal4VP16/UAS-ECFP double-transgenic mice[END_REF]), Csf1r-mApple (MacApple, [START_REF] Hawley | Csf1r -mApple Transgene Expression and Ligand Binding In Vivo Reveal Dynamics of CSF1R Expression within the Mononuclear Phagocyte System[END_REF]) and PyMT-ChOVA mice were intercrossed to generate MacBlue x Cx3cr1 EGFP/+ x PyMT-ChOVA and MacBlue x Cx3cr1 EGFP/+ x MacApple x PyMT-ChOVA mouse strains. These strains were bred at the Centre d'Exploration Fonctionnelle Pitié-Salpêtrière animal facility. To monitor TAM accumulation over time, between 14 and 40-week old mice were used. All mice were maintained under SPF conditions at 22°C and used after breast tumor development. OTI-1 RAG-2 -/-mice were kindly provided by Sebastian Amigorena (Institut Curie, France). All experiment protocols were approved by the French animal experimentation and ethics committee and validated by Service Protection et Santé

Quantification and Statistical Analysis

All statistical analyses were performed with Graphpad Prism 7 for flow cytometry and imaging, and R for scRNAseq analysis. Multigroup analysis of variances were performed, and one-way or two-way ANOVA tests were performed followed by Turkey post test for Gaussian distribution or Kruskall-Wallis followed by Dunn's multiple comparisons. For simple comparison analysis, Mann-Whitney was performed to compare nonparametric distribution. For correlation analysis, Pearson correlation was calculated *, P < 0.05; **, P < 0.01; ***, P < 0.001. Sample sizes are indicated in each figure legend.

SUPPLEMENTARY MOVIE

Click here to access/download;Figure;