Transcriptional Profiles and Stromal Changes Reveal Bone Marrow Adaptation to Early Breast Cancer in Association with Deregulated Circulating microRNAs Free

Carcinogenesis is a complex and dynamical process centered on the transforming cells and in need of surrounding stromal and immune microenvironment (1, 2). Besides local tissue modifications, cancer progression induces a profound and systemic reshaping of the host immune functions toward an immunosuppressive state. This systemic influence has been mostly investigated in transplantable tumor models, in which the expansion of bone marrow–derived myeloid cells endowed with immunosuppressive activity has been associated with cancer size (3–6) but also with correlative evidence in patients with cancer (7). Despite such consolidated view of cancer as a systemic disease, little is known about the capacity of nascent tumors to interlace a cross-talk with the bone marrow toward the instruction of a tumor-educated hematopoiesis and even less known is how early in transformation the bone marrow can sense and respond to a distant tumor. Indeed, despite very elegant single-cell resolution studies mapping the bone marrow transcriptional landscape at homeostasis and stress-induced haematopoiesis, such as in case of chemotherapy (8), a transcriptional profile of the bone marrow hematopoietic compartment along carcinogenesis is still missing.

To fill this gap, we investigated, by gene expression profile (GEP), the modifications occurring in the bone marrow hematopoietic environment along discrete stages of tumor development, that is, from premalignant to invasive stages, in the spontaneous MMTV-NeuT mouse model of mammary carcinogenesis (9). The finding of a transcriptional reprogramming of the hematopoietic compartment already at preinvasive stage of neoplastic progression has been correlated with modifications in both the cellular composition and the stromal architecture of the bone marrow. Functional study has been performed testing the top upmodulated gene in the bone marrow at early disease stage for the capacity of reproducing features of cancer-associated myelopoiesis potentially contributing to tumor growth.

Moreover, we tested the hypothesis that circulating molecules could message the cross-communication between nascent cancers and bone marrow. To this aim, we profiled plasma miRNAs along disease progression because of their regulatory role in many biological programs (10). Changes in miRNAs expression during neoplastic progression, if correlated with bone marrow transcriptional changes, could represent relevant noninvasive biomarkers of disease.

Our study offers a new insight into bone marrow hematopoiesis adaptation to the incipient peripheral transformation and starts defining some molecular traits of this dynamic process.

All animal studies were approved by Institutional Committee for Animal Welfare and by the Italian Ministry of Health and performed in accordance with National Law D.lgs 26/2014 (authorization INT_05_2012). BALB-NeuT female mice expressing the transforming rat oncogene c-erbB2 (Her-2/neu) under the mouse mammary tumor virus (MMTV) promoter spontaneously develop mammary carcinoma with a well-defined tumor progression (9, 11). PyMT mice (C57BL/6 genetic background) express the oncogene polyoma middle T antigen under the MMTV promoter (12). HuHer2Δ16 mice (FVB genetic background), kindly provided by S. Pupa (Fondazione IRCCS Istituto Nazionale dei Tumori, Milan, Italy), carry, under MMTV promoter, the human HER2 splice variant lacking exon-16 (13).

Bone marrow cells were obtained by flushing femurs and tibiae with medium using a 22-gauge needle attached to a 1-mL syringe. Five hundred microliters of blood was collected by intracardiac puncture with 20 μL of EDTA 0.5 mol/L pH 8 for FACS analysis and for plasma separation within 2 hours.

Bones and mammary glands were fixed with 10% neutral-buffered formalin overnight, washed in water, bones decalcified, and paraffin-embedded. Four-micrometer-thick sections were used for both histopathology and immunostaining. IHC was performed using the horseradish peroxidase method as reported previously (14). See Supplementary Materials and Methods for a list of antibodies.

Confocal microscopy analysis for CD29, SPARC, and c-Kit has been performed on frozen, OCT-embedded bone marrow samples. Sections were fixed and incubated with mAbs to SPARC (AF952, R&D Systems) followed by its secondary anti-goat Alexa 546-conjugated Ab (Invitrogen). CD29-APC conjugated (eBioHMb1-1; eBioscience) and c-Kit FITC-conjugated (2B8; eBioscience) Abs have been added as third step after a block in a 10% FCS solution. Images were acquired under a Leica DM4 B optical microscope equipped with a Leica DFC450 digital camera.

Single-cell suspension preparations were treated with ammonium chloride potassium lysis buffer to remove red blood cells, washed, and incubated with specific antibodies (see Supplementary Materials and Methods for a list of antibodies). Samples were acquired using a BD LSR II Fortessa instrument (Becton Dickinson) and analyzed with FlowJo software (version 8.8.7, TreeStar). All samples were analyzed in single.

To construct tissue-specific lentivectors, we modified the self-inactivating lentiviral vector pRRL.GFP.WPRE.hPGK.NGFR.GATA-1 (backbone pRRLCMVGFPsin-18; a kind gift from Dr. G. Ferrari, IRCCS San Raffaele Scientific Institute, Milan, Italy; ref. 15) by replacing the GATA-1 promoter with myeloid-specific hCD68sh (16) and then the GFP sequence with Atf3 cDNA. A third-generation packaging system was used to produce viral particles. Lentiviral stocks were produced in 293T cells by Ca₃PO₄ cotransfection of the plasmids as described previously (17). Bone marrow lineage-negative (Lin⁻) cells were purified using the Mouse Lineage Cell Depletion Kit (Miltenyi Biotec), following the manufacturer's instructions. Lin⁻ cells isolated from donor mice were prestimulated overnight in StemSpan SFEM (STEMCELL Technologies) supplemented with 100 ng/mL SCF, 100 ng/mL Flt3L, 50 ng/mL TPO, 20 ng/mL IL3 and then lentiviral particles added at a multiplicity of infection (MOI) of 100. After 24 hours, cells were washed, checked for efficiency of transduction by flow cytometry and i.v. injected, at the dose of 5×10⁁5 cells/mouse, into lethally irradiated (split dose 500+500 RAD) mice.

Total RNA was extracted from bone marrow cells using TRIzol reagent (Life Technologies) following the manufacturer's instruction, treated with DNase 1 (Qiagen), and quantified by spectrophotometry (ND-2000c; NanoDrop Products). RNA integrity was verified using the RNA 6000 Nano Kit (Agilent Technologies). For GEP, RNA was reverse transcribed, labeled with biotin, and amplified overnight using the TotalPrep RNA Amplification Kit (Illumina). Biotinylated complementary RNA was hybridized to mouseWG-6 v2.0 Expression BeadChips (Illumina). Arrays were scanned with Illumina BeadArray Reader, and raw data obtained using Illumina BeadStudio v3.3.8 were processed using the lumi package (18) from R/Bioconductor (19). Raw data were log₂ transformed, normalized using robust spline normalization, and filtered: Only probes with a detection P < 0.01 in at least one sample were selected for downstream analyses. Multiple probes representing the same gene were collapsed selecting the probe with the highest detection rate across samples. In case of ties, the probe with the highest interquartile range was selected. Differentially expressed genes (DEG) were identified using the limma package (20). P values were adjusted for multiple tests using the Benjamini–Hocheberg FDR. Genes with a fold change ≥1.5 and FDR < 0.05 were considered statistically significant. Hierarchical clustering was applied with average linkage and 1-Pearson correlation coefficient as distance metric.

Expression data were deposited in the NCBI Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/geo/) with the accession number GSE117071.

Gene set enrichment analysis (GSEA; ref. 21) was carried out in preranked mode using the c2_all collection from MSigDB database (http://software.broadinstitute.org/gsea/msigdb/index.jsp). Gene sets with an FDR < 0.05 were considered statistically significant. Functional overrepresentation analysis of DEGs was carried using the topGO package with Gene Ontology (GO) biological process terms, and QIAGEN's Ingenuity Pathway Analysis (QIAGEN, www.qiagen.com/ingenuity). Single sample GSEA (ssGSEA; ref. 22) was performed using Kyoto Encyclopedia of Genes and Genomes (KEGG) gene sets pathways. Limma analysis was done between BALB/c and NeuT mice at each time point and the comparisons were considered significant with FDR < 0.05.

To estimate the relative amount of immune cells in the bone marrow samples, we applied ImmQuant, a computational algorithm that allows to infer the relative composition of different immune cell types in mouse tissues according to mRNA expression data (23). We have run ImmQuant on BALB/c and NeuT at 12 and 24 weeks, selecting BALB/c as the control group. Significant modulations of immune cell-type quantities were calculated through two-sample t test (P < 0.05).

Total RNA was extracted from plasma samples using the column-based system miRNeasy Mini Kit (Qiagen), slightly modified by Exiqon as described before (24).

To assess the expression levels of miRNAs Agilent SurePrint G3 Mouse miRNA 8 × 60K arrays designed on miRBase 19.0 were used as described previously (24). Microarrays were scanned with an Agilent SureScan, and raw data were collected using Agilent's Feature Extraction software v10.7. Raw data were background-corrected and log₂ transformed using the rmaMicroRna function of the AgiMicroRna package. Because of the presence of a batch effect, miRNA expression data were corrected using ComBat (25). miRNAs detected in at least 50% of samples, as indicated by the gIsGeneDetected parameter given by the Feature Extraction software, were kept for further analyses. miRNAs differentially expressed at each time point were identified using the limma package. Multiple-testing correction was performed using the Benjamini–Hochberg FDR and miRNAs with a nominal P < 0.05 were considered significant. Expression data were deposited in the GEO with the accession number GSE117071.

RNA (1 μg) was reverse transcribed with the High Capacity cDNA Reverse Transcription Kit from Applied Biosystems according to the manufacturer's instructions. Real-time PCR (RT-PCR) reactions were performed on ABI Prism 7900 HT (AB) using the TaqMan Fast Advanced PCR MasterMix (Applied Biosystems) according to the manufacturer's instructions. Gene expression levels were normalized through the comparison with mouse Gapdh expression (see Supplementary Materials and Methods for a list of RT-PCR primers). Results were obtained using the comparative C_t method (26).

For validation of array data, plasma RNA was reverse transcribed with the TaqMan miRNA Reverse Transcription Kit and miRNA specific stem-loop primers (Applied Biosystems) in a small-scale RT reaction as described in ref. 27. The RT-PCRs were conducted on an ABI Prism 7900 HT using the TaqMan Fast Advanced Master Mix (Applied Biosystems). Data were analyzed with SDS software version 2.2.2 (Applied Biosystems) with automatic cycle-threshold (C_t) settings for assigning baseline and threshold. Synthetic C. elegans miRNA, cel-miR-39 (IDT), was added after sample denaturation and used as technical control of sample processing. The new cohort of plasma samples was analyzed with custom RT and preamplification pools using TaqMan miRNA assays according to manufacturer's instructions.

To assess whether genes downmodulated in the bone marrow of NeuT mice are targeted by differentially expressed miRNAs, we used miRWalk2.0 selecting 5 algorithms (miRWalk, PITA, RNA22, miRanda, and Targetscan). Among the putative target genes, we considered only those predicted by at least 3 algorithms. Such gene lists were submitted to the DAVID (Database for Annotation, Visualization and Integrated Discovery) using the Functional Annotation Tool to discover whether those genes belong to known pathways.

The BLASTN algorithm implemented in miRBase database (28) was used to search for predicted human homologous sequences of the differentially expressed miRNAs, excluding those alignments with an e-value ≥ 0.05. GEO dataset GSE83270, which includes miRNA profiles from the blood of 6 patients with breast cancer and 6 healthy controls, was considered for comparative analysis. Preprocessed data were imported in R software and differentially expressed miRNAs between patients with cancer and healthy donors were identified by two-sided Wilcoxon test with P < 0.05.

Microarray miRNA profiles of human breast tumors included in the METABRIC cohort (29) were obtained upon request from the European Genome-phenome Archive, with accession number EGAD00010000438. Clinical data for METABRIC patients were retrieved from cBioPortal (http://www.cbioportal.org). Normalized miRNA expression data were reannotated to miRBase version 20. Luminal-A and luminal-B patients according to the PAM50 classification were selected for survival analysis. For each miRNA, samples were classified as high- or low-expressing according to the median.

Data analysis was performed using Prism software (GraphPad Software, Inc.). Results are expressed as means ± SD, as specified. Statistical analysis of continuous variables was performed using a two-tailed Student t test with confidence intervals of 95%. Data were considered significantly different at P < 0.05 (*, P < 0.05; **, P < 0.01; ***, P < 0.005). For survival analysis in METABRIC cohort, Kaplan–Meier curves, log-rank test, and Cox proportional hazards model were used, as implemented in the survival R package. For multivariate Cox model in addition to miRNA expression, we included age at diagnosis, tumor size, tumor grade, number of positive lymph nodes, and PAM50 subtype.

Using GEP, we interrogated the spontaneous mammary adenocarcinoma MMTV-NeuT (hereafter NeuT) model for the modifications occurring in the bone marrow transcriptional landscape along discrete disease stages of mammary transformation.

We profiled bone marrow cells of NeuT mice at 6, 12, and 24 weeks of age, roughly representative of premalignant, preinvasive, and invasive stages of transformation in comparison with age-matched normal siblings (Fig. 1A). Premalignant and preinvasive stages of mammary gland transformation were scored by histopathologic analysis performed on the primary lesions of each mouse (see Supplementary Table S1 and Supplementary Materials and Methods).

Most of the DEGs were from the paired comparison at 24 weeks, with 141 genes upmodulated (fold change ≥ 1.5, FDR ≤ 0.05) and 233 genes downmodulated (fold change ≤ −1.5, FDR ≤ 0.05; Fig. 1B and C; Supplementary Table S2). At the preinvasive stage, that is, 12 weeks, only 10 genes were significantly upmodulated in NeuT mice (Fig. 1B and D). At 6 weeks, the premalignant stage, only one gene, Peg3 (paternally expressed 3), was significantly upmodulated in transgenic mice. As Peg3 resulted strikingly upmodulated in the bone marrow of NeuT mice at every time point, we tested whether this upregulation could be a consequence of the transgene integration. TLA (targeted locus amplification) technology demonstrated that the NeuT-expressing vector had integrated within a region of chromosome 7 where Peg3 is located, likely causing a genomic rearrangement. In light of this evidence and of other investigations (Supplementary Materials and Methods; Supplementary Fig. S1), Peg3 was not included in the subsequent analyses.

In the bone marrow of 12-week-old NeuT mice, the top upregulated DEGs were related to innate immune response and included Atf3 (activating transcription factor 3), Trem1 (triggering receptor expressed on myeloid cells 1), Tnfaip3 (TNF, alpha-induced protein 3), and Ncf2 (neutrophil cytosolic factor 2). The upregulation of these genes was technically validated by RT-PCR (Supplementary Fig. S1D).

The biological networks involved in the transcriptional reprogramming were assessed by GO on the genes differentially modulated at 24 weeks, as their number was sufficient for such analysis. Most of the downmodulated genes in NeuT mice were related with immune responses, T-cell differentiation and activation and B-cell receptor signaling pathway, while upregulated genes were mostly related with innate immune response, leukocyte activation and response to danger signals (Fig. 1E; Supplementary Table S2). SsGSEA (22) confirmed a clear trend of changes in immune-related networks during disease progression with a downmodulation of adaptive immune response programs, such as T-cell receptor signaling and antigen presentation pathways, and an upregulation of those relative to immunosuppression and inflammation (TGFβ and TNF networks; Supplementary Fig. S1E).

To focus on the bone marrow changes more closely associated with disease progression, we compared the NeuT transcriptional profiles longitudinally (6, 12, and 24 weeks; Fig. 2A), filtering out DEGs in the bone marrow of their control siblings of the same age. This approach should exclude the genes related to age and/or environmental conditions. The longitudinal analysis indicated a number of significantly modulated genes, of which, 40 up- and 33 downregulated already in the 12 versus 6 weeks comparison (Fig. 2B; Supplementary Table S3). GO analysis indicated enrichment in terms related to immune responses and extracellular matrix (ECM) organization (Fig. 2C). In line, Sparc (secreted protein acidic and rich in cysteine), a matricellular protein involved in collagen deposition and ECM organization, emerged among the most downmodulated genes (FC, −2.56; FDR, 0.0111), and its downmodulation was confirmed at protein level by IHC (Fig. 2D and E). Double immunofluorescence (IF) analysis showed that SPARC expression in the bone marrow mostly colocalizes with CD29⁺ stromal cells and megakaryocytes (Fig. 2F; Supplementary Fig. S2) according to a previous report (14).

DEGs between NeuT mice at late disease stage and control siblings (with a fold change ≥ 2) were used to derive a cancer-adapted hematopoiesis signature that was tested for its ability to correctly distinguish NeuT mice from age-matched controls, according to their bone marrow transcriptional profile, before tumor was manifested. Strikingly, on unsupervised hierarchical clustering, this signature neatly separated NeuT from control samples (Fig. 2G). To further validate this signature independently from the driving oncogene, we profiled bone marrow cells from two additional mouse mammary tumor models, namely MMTV-PyMT and MMTV-huHer2Δ16, at the infiltrative cancer stage, and age-matched syngeneic controls, C57BL/6 and FVB, respectively. GSEA showed that genes included in the signature were coherently enriched among the genes up- and downmodulated in the two additional transgenic models (Fig. 2H), suggesting that the transcriptional reprogramming of the bone marrow hematopoiesis can be a common feature of transformation toward breast cancer, independently from the mouse background and from the driving oncogene.

To investigate which immune cell types were mainly associated with bone marrow transcriptional reprogramming, we applied to our GEP data the computational algorithm ImmQuant, which allows to infer the relative composition of different immune cell types in mouse tissues according to mRNA expression data (23). The deconvolution analysis indicated a significant increase in neutrophils and decrease in B cells at both 12- and 24-week time points. In addition, at late-stage disease, ImmQuant predicted the enrichment of classical monocytes and macrophages. B-cell decrease involved different B-cell subpopulations, including precursor cells (such as preB and proB). At early stages, in addition to Ly6G⁺ neutrophils, CD103⁺ dendritic cells and CD25⁺ T regulatory cells (Treg) were also increased (Fig. 3A and B).

To test the actual correspondence between data inferred from the bone marrow GEP and changes in the bone marrow hematopoietic composition, we performed histopathologic and flow cytometry analysis along disease progression beginning with the overt infiltrative stage. Bone marrow parenchyma showed a marked granulocytic myeloid hyperplasia, characterized by an enrichment in morphologically immature granulocytic elements, associated with a contraction of erythroid precursors and lymphoid elements (Supplementary Fig. S3A), all features confirmed by IHC for MPO (myeloid cells), Ter-119 (erythroid cells), and Pax-5 (B cells; Supplementary Fig. S3B and S3C) and by flow cytometry analysis (Supplementary Fig. S3D and S3F). In line with the expansion of the myeloid compartment, we found a significant increase in the GMP (granulocytic–monocytic progenitors) fraction of Lin⁻ c-kit⁺ hematopoietic precursors (Supplementary Fig. S3E and S3F).

To define how early in transformation these changes occur, bone marrow at earlier stages has been analyzed. At 12 weeks, a clear-cut expansion of granulocytic myeloid cells and a reduction of erythroid and lymphoid elements were detected in the bone marrow of NeuT mice. The same modifications, although less prominent, were also detected in the marrow of 6 weeks NeuT mice (Fig. 3C–F), in which the mammary glands showed moderate-to-severe epithelial dysplasia.

In light of the emerging role of the different T-cell subsets in the bone marrow in advanced/metastatic settings (30), we tested whether changes in this compartment could be detected at early disease stage (12 weeks). Multiparametric flow cytometry showed an overall decrease of both CD4 and CD8 T cells, with an increased percentage of Foxp3⁺ regulatory T cells (Fig. 3G). In addition, both CD4 and CD8 T cells shifted toward the effector memory subset at the expense of naïve and central memory cells (Fig. 3H).

Overall, these data indicated that the transcriptional reprogramming of the bone marrow detected already at early stage of carcinogenesis is associated with significant changes in its cellular composition.

We next investigated whether the changes in the bone marrow cellular composition were associated with modification of the stromal architecture beginning at early stage of transformation.

Histopathologic analysis showed alterations in the stromal architecture of the bone marrow at both 12- and 6-week time points. The density and localization of Nestin⁺ mesenchymal stromal elements and the expression of the prototypical stromal-derived chemotactic factor CXCL12 were both increased in the hematopoietic interstitium of transgenic mice than control siblings (Fig. 4A and B). Indeed, while in the hematopoietic lacunae of control mice Nestin⁺ and CXCL12-expressing elements were mostly confined to the vascular area, in NeuT mice the expression of Nestin and CXCL12 were more prominent in the perivascular and interstitial areas. Such modifications were quantified by an ad hoc software analysis applied to immunostained bone marrow sections (Supplementary Fig. S4A). The relocalization of Nestin⁺ meshwork and CXCL12 expression suggested an alteration in the hematopoietic niche and was associated with a change in the expression and immunolocalization of the CXCL12 receptor CXCR4, which was consistently increased in myeloid cells of 6- and 12-week NeuT mice within interstitial areas (Supplementary Fig. S4B). The actual cooccurrence of CXCL12/CXCR4 induction in the same hematopoietic environment was confirmed by double marker immunofluorescence analysis on 12-week bone marrow samples showing the coherent modulation of CXCL12 and CXCR4 in the hematopoietic interstitium (Supplementary Fig. S4C).

Further supporting the association between early-occurring NeuT stromal changes and modifications in the topology of hematopoietic populations, confocal microscopy analysis showed a downmodulation of SPARC expression in mesenchymal CD29⁺ elements lining the endosteal niche, along with an increase and relocalization of cKit⁺ precursor cells from osteoblastic niche to interstitial areas (Fig. 4C).

These results indicate that the hematopoietic switch observed in the bone marrow hematopoiesis along tumor progression is sustained by detectable changes in the stromal niche organization.

The transcription factor Atf3, a major molecular hub of inflammatory response to stress and danger sensing (31), was the top upmodulated gene in the comparison between NeuT and control sibling bone marrow of 12 weeks of age (Supplementary Table S2; Fig. 1D). Its expression in the bone marrow perfectly correlates with the pathologic disease score of mammary lesions in individual mice (Fig. 5A). Transcriptional upmodulation of Atf3 correlated with protein expression (Fig. 5B). Notably, in the bone marrow of NeuT mice, Atf3 showed a clear nuclear localization, suggestive of its activation, whereas in the bone marrow of control siblings, it is mainly cytoplasmic (Fig. 5C). To confirm that the observed upmodulation of Atf3 expression was not only due to the expansion of the myeloid compartment, we checked the ratio between Atf3 and CD11b expression, calculated for every single mouse of the GEP cohort and confirmed that Atf3 upmodulation was not only due to the expansion of the myeloid compartment, but also to increased expression at single cell level (Fig. 5D, left). In addition, myeloid cell subtype purification by microbeads indicated that Atf3 overexpression in NeuT bone marrow was mainly in the Gr1^dimLy6G⁻ subset (Fig. 5D, right).

Mechanistically, we tested whether the upmodulation of Atf3 in the bone marrow of NeuT mice could be functionally linked with the expansion of the myelopoietic compartment. Lin⁻ bone marrow cells from naïve BALB/c mice were transduced with a lentiviral vector expressing Atf3 under the myeloid-specific CD68 promoter and transplanted into lethally irradiated NeuT and control siblings. After 10 weeks, mice were sacrificed and the bone marrow collected for histopathology, IHC, and flow cytometry. Although the efficiency of Lin⁻ transduction was suboptimal, ranging between 15% and 30% of total bone marrow cells, mice with the highest transduction level showed the most relevant expansion of the myeloid compartment (MPO⁺ cells), paralleled by a contraction of erythroid and B-cell pools (Fig. 5E–I). Notably, bone marrow of WT mice receiving Atf3-transduced Lin⁻ cells showed enhanced hematopoiesis similar to that of nontransplanted NeuT mice. This result suggests a general effect of Atf3 in unleashing myelopoiesis, regardless of the presence of the NeuT oncogene in recipient mice. The same Atf3-transduced Lin⁻ cells transplanted into NeuT recipients further exacerbated oncogene-induced myelopoiesis. We did not observe significant modifications in Nestin/CXCL12 and CXCR4 expression (Supplementary Fig. S4D) despite the "enforced" myelopoiesis. This suggests that Atf3 overexpression bypasses the need of stromal niche modification for the expansion of myeloid cells.

Of note, histopathologic analysis of the mammary lesions 10 weeks after bone marrow transplantation showed a higher tumor burden in NeuT mice receiving Atf3-expressing than control Lin⁻ cells (Fig. 5J). These data point to Atf3 upregulation in the myeloid compartment as a potential positive regulator of tumor-supportive central hematopoiesis.

Having demonstrated that stromal and hematopoietic elements of the bone marrow parenchyma act as early sensor of malignant mammary transformation, we looked at miRNAs as soluble mediators potentially responsible for bone marrow transcriptional alterations. Profiling of plasma miRNAs revealed that 80 miRNAs were differentially expressed between NeuT and age-matched siblings at the overt infiltrative cancer stage (24 weeks), while 28 miRNAs were already differentially upmodulated at the preinvasive stage (12 weeks; Fig. 6A; Supplementary Table S4). Among the latter group, we identified miRNAs belonging to the miR-29 (miR-29a-3p; miR-29c-3p) and miR-23 (miR-23a-3p; miR-23b-3p) clusters, which are known to modulate the extracellular matrix synthesis (32) and the CXCL12/CXCR4 axis (33), respectively.

miRNAs differentially expressed at the preinvasive stage were validated by RT-PCR on the same plasma used for microarrays as well as in an independent cohort of NeuT mice and age-matched control siblings (Supplementary Fig. S5A; Fig. 6B). Notably, the expression of miR-29a and miR-23a correlated with the tumor burden in the mammary glands from 12-week NeuT mice (Fig. 6C), fueling the hypothesis of a preferential cancer influence over circulating miRNA levels.

To further investigate the origin of the above miRNAs, we have analyzed their expression in epithelial cells from normal mammary glands and from established tumor nodules, as well as in the interstitial fluid of mammary glands from 12 weeks old NeuT or control siblings.

miR-23b-3p, -27a-3p, and -671-5p are more expressed in neoplastic cells than normal mammary epithelium, and more present in the interstitial fluid from early mammary glands of NeuT than control mice (Fig. 7A). Other miRNAs, such as miR-29a-3p, -29b-3p, and -29c-3p, are not differentially expressed between normal and neoplastic mammary cells, neither in the interstitial fluid from NeuT or control mammary glands, suggesting their release from circulating immune cells.

We then assessed whether differentially expressed circulating miRNAs could target directly genes downregulated in the bone marrow of NeuT mice. Using miRWalk2.0 with 5 algorithms (miRWalk, PITA, RNA22, miRanda, Targetscan), we identified 14,431 predicted gene targets (at least by 3/5 algorithms) of any of the 28 miRNAs differentially expressed at 12 weeks. Of these genes, 171 were significantly downmodulated in the bone marrow of NeuT mice and were targeted, on average, by 6 miRNAs. Similarly, we identified 16,970 predicted targets of the 80 miRNAs differentially expressed at 24 weeks, of which, 209, targeted on average by 13 miRNAs, were significantly downregulated in the NeuT bone marrow (Supplementary Fig. S5B). KEGG pathway analysis revealed that most of the 209 predicted targets belong to B-cell receptor signaling and antigen processing and presentation pathways (Fig. 7B), in line with the alterations observed in bone marrow hematopoiesis. These results suggest that circulating miRNA imbalance, associated with breast cancer development, can be detected even in preinvasive conditions, offering correlates with the bone marrow transcriptional and phenotypic switch of cancer-adapted hematopoiesis.

To investigate the translational significance of these findings, we compared our list of circulating miRNA with the profiles of patients with breast cancer versus healthy controls in one of the very few publicly available datasets (GSE83270; ref. 34). Among the human–mouse homologous miRNAs (39/94), 11 were significantly upregulated in the blood of patients with breast cancer compared with healthy controls (Fig. 7C).

Moreover, we assessed whether the expression of any of these miRNAs identified in the NeuT model was associated with overall survival in patients with luminal breast cancer included in the METABRIC dataset. According to univariate Cox model, we found 5 miRNAs significantly associated (FDR < 0.05) with worse and 5 with better disease-specific survival (Supplementary Fig. S5C). In a multivariate Cox model (including age at diagnosis, tumor size, tumor grade, number of positive lymph nodes, and PAM50 subtype as covariates), 3 miRNAs (hsa-miR-202-3p, hsa-miR-29a-3p, and hsa-miR-125a-3p) retained their significance, indicating that the prognostic relevance of these miRNAs is independent of the other clinical–pathologic and molecular parameters tested.

Bone marrow stromal cells act as sensors of immunologic stress signals from peripheral tissue. As for infections and autoimmunity, the presence of solid tumors induces a chronic state of inflammation that involves the bone marrow. Among bone marrow responses to a distant growing tumor, the expansion and release of myeloid cells with suppressive activities is certainly the most studied (35), in both mouse models (3–5, 36) and patients with cancer (7). Although this phenomenon has been extensively described in case of advanced diseases, the question on how early from initial transformation the bone marrow can sense and react to a nascent tumor remains largely unexplored. Here, we provide evidence of significant changes occurring in the bone marrow already at early disease stage, when only preinvasive lesions are detectable in the mammary glands of the transgenic NeuT model of breast cancer. Bone marrow alterations, which become more evident with disease progression, include a transcriptional reprogramming of the hematopoietic compartment, associated with changes in specific immune cell subsets, and modifications in the stromal architecture.

Among the biological programs characterizing the bone marrow transcriptional signature at invasive disease stage we found innate immune responses and defense responses to external stimuli as upmodulated, whereas B-cell and T-cell lymphoid activation/differentiation programs were coherently downregulated. These same biological programs were found enriched in the bone marrow transcriptome of NeuT mice with preinvasive disease, suggesting that bone marrow hematopoietic homeostasis is perturbed early during carcinogenesis. Supporting this hypothesis, the application of this cancer-adapted hematopoiesis signature allowed unsupervised separation of transgenic mice from their sibling wild-type counterparts at very early disease stages (6 and 12 weeks).

Deconvolution analysis of the bone marrow gene expression profiles with ImmQuant tools allowed to infer alterations in specific immune cell subsets, in particular the expansion of myeloid subsets (mainly granulocytes), at the expense of the B-cell compartment, as well as an increase in Tregs. Confirmation of the in silico evidence came from flow cytometry and histopathologic analysis, the latter also highlighting a contraction of the Ter119⁺ erythroid fraction. The alteration of the myeloid-to-erythroid ratio of the bone marrow hematopoietic parenchyma has previously been reported in the spontaneous model of invasive breast cancer (PyMT), associated with the expansion of Ly6G⁺ cells (5) and in transplantable mammary tumor cell lines in which overt tumor growth was associated with increased myeloproliferation, impaired erythropoiesis, and loss of early bone marrow progenitor cells (37). These common traits have been anticipated to the preinvasive disease stage in our study, which offers new clues on the earliest timing in which the bone marrow senses and starts responding to a distant transformation. Flow cytometry analysis of the bone marrow T-cell compartment also confirmed the predicted enrichment in Tregs at early disease stage that is consistent with the enhanced expression of CXCL12, accordingly to the relevance of CXCR4/CXCL12 signaling for Tregs trafficking in the bone marrow (30, 38).

The functional relevance of genes upmodulated in the bone marrow of NeuT mice at the preinvasive stage comes from the analysis of Atf3, the top upregulated gene. Atf3 is an immediate early-response gene upregulated upon sensing of stress signals by myeloid cells, for regulation of cell–cell communication and migration programs and for regulation of the inflammatory response via NF-κB (39, 40). Myeloid-specific overexpression by gene transduction into Lin⁻ precursor and transplantation into BALB/c mice (void of tumor) demonstrates that ATF3 forced expression phenocopies many aspects of tumor-induced bone marrow myelopoiesis, whereas the same transplant performed into NeuT mice showed exacerbation of mammary tumor progression. ATF3 has been already implicated in breast cancer progression and metastasis through its role in a TGFβ feed-forward loop in malignant epithelial cells (31). Its involvement in cancer-supportive myeloid populations has been demonstrated in the PyMT orthotopic model in which LysM-Cre mediated deficiency of ATF3 altered tumor-associated macrophages toward a proinflammatory phenotype (41). These pieces of evidence claim for further investigation on the role of ATF3 in the functional reprogramming of myeloid cells and in early hematopoietic adaptation to cancer. Nevertheless, our data suggest that in bone marrow Atf3 regulation in the myeloid compartment might be at the level of stromal niches and that its enforced expression bypasses such regulation.

Bone marrow stromal niches are essential for the maintenance and regulation of hematopoietic stem cells and of their lineage differentiation (42, 43). Also, stromal remodeling occurs in response to stress conditions or in case of hematologic disorders (8, 14, 44). Much less is known on bone marrow stromal niche alterations associated to distant solid tumor, independently from bone marrow metastatic colonization. Recently, a cross-talk between lung tumors and bone marrow osteoblasts has been demonstrated, which is functional to the expansion of SiglecF⁺ neutrophils with a role in tumor growth, immune escape, and metastasis (36). In the bone marrow of NeuT transgenic mice, changes in the topographical expression of CXCL12 stromal chemokine and in the overall distribution of Nestin⁺ mesenchymal cells begin at preinvasive stage. Modulation of CXCL12 in the interstitial stromal meshwork was associated with the accumulation of CXCR4-expressing myeloid elements in the bone marrow of transgenic mice. The modulation of CXCL12/CXCR4 gradients has been reported to play a central role in hematopoietic adaptation to peripheral stimuli (45) and in the attraction of suppressive myeloid cells associated with cancer development (46). CXCL12 modulation in the bone marrow stroma along cancer progression has not been investigated so far. Such alteration in the CXCL12/CXCR4 axis could represent one mechanism involved in the modified bone marrow hematopoietic composition of NeuT mice.

Associated to alteration of the CXCL12/CXCR4 axis, the bone marrow of NeuT mice shows downmodulation of the matricellular protein SPARC, which decreases along disease progression. The lower amount of SPARC in CD29⁺ mesenchymal elements lining the endosteal niche likely contributes to the spatial relocalization of cKit⁺ precursors to interstitial areas of the bone marrow parenchyma. SPARC activity is epigenetically regulated by the miR-29 family (32), which includes miR-29a-3p and miR-29c-3p that we have found overexpressed in the plasma of NeuT mice. Besides SPARC, the predicted target list of the miRNAs significantly induced in the blood of NeuT mice includes genes involved in B-cell receptor signaling and function, accordingly found downmodulated in their bone marrow. This evidence represents a potential link between peripheral transformation and bone marrow hematopoietic adaptation through circulating miRNAs. Although the strength of this connection needs further investigation, the finding of circulating miRNA deregulation already at the preinvasive stage of breast cancer progression could represent an important hint for translation in the human setting. Their potential translational relevance is also suggested by the differential expression of 11 human homologs of our circulating miRNAs in a publicly available dataset of patients with breast cancer compared with healthy controls and the prognostic performance of 5 of them in the METABRIC human luminal breast cancer cohort.

Altogether our findings offer a new perspective on the systemic rewiring of the immune system during carcinogenesis, defining bone marrow hematopoietic adaptation as a progressive event preceding overt malignancy and articulated around a core of innate/inflammatory response programs.

No potential conflicts of interest were disclosed.

Conception and design: C. Chiodoni, T.A. Renzi, C. Tripodo, M.P. Colombo

Development of methodology: M. Perrone, M.P. Colombo

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): C. Chiodoni, T.A. Renzi, M. Perrone, A.M. Tomirotti, S. Sangaletti, L. Botti, C. Tripodo

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): C. Chiodoni, V. Cancila, T.A. Renzi, A.M. Tomirotti, M. Dugo, M. Milani, M. Marrale, C. Tripodo

Writing, review, and/or revision of the manuscript: C. Chiodoni, T.A. Renzi, M. Milani, C. Tripodo, M.P. Colombo

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): V. Cancila, L. Bongiovanni

Study supervision: C. Chiodoni, C. Tripodo, M.P. Colombo

The research leading to these results has received funding from AIRC under 5 × 1000 Special Program "Tumor microenvironment related changes as new tools for early detection and assessment of high risk disease" #12162 (M.P. Colombo and C. Tripodo as group leaders) and from Italian Ministry of Health (Ricerca Corrente IZS SA 04/15). T.A. Renzi was supported by AIRC fellowship "Terme di Sirmione" # 19425. M. Perrone is a student registered with the Open University (Milton Keynes, United Kingdom). We thank the Platform of Integrated Biology for microarray experiments, the Unit of Bioinformatics and Biostatistics, the Animal Facility at Fondazione IRCCS Istituto Nazionale dei Tumori, and Ester Grande for administrative support. This study makes use of data generated by the Molecular Taxonomy of Breast Cancer International Consortium funded by Cancer Research UK and the British Columbia Cancer Agency Branch.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.