Evolution of fibroblasts in the lung metastatic microenvironment is driven by stage-specific transcriptional plasticity

Mortality from breast cancer is almost exclusively a result of tumor metastasis, and lungs are one of the main metastatic sites. Cancer-associated fibroblasts are prominent players in the microenvironment of breast cancer. However, their role in the metastatic niche is largely unknown. In this study, we profiled the transcriptional co-evolution of lung fibroblasts isolated from transgenic mice at defined stage-specific time points of metastases formation. Employing multiple knowledge-based platforms of data analysis provided powerful insights on functional and temporal regulation of the transcriptome of fibroblasts. We demonstrate that fibroblasts in lung metastases are transcriptionally dynamic and plastic, and reveal stage-specific gene signatures that imply functional tasks, including extracellular matrix remodeling, stress response, and shaping the inflammatory microenvironment. Furthermore, we identified Myc as a central regulator of fibroblast rewiring and found that stromal upregulation of Myc transcriptional networks is associated with disease progression in human breast cancer.


50
Breast cancer continues to be one of the leading causes of cancer related death in women, 51 and mortality is almost exclusively a result of tumor metastasis. Advanced metastatic 52 cancers are mostly incurable and available therapies generally prolong life to a limited extent.

53
It is increasingly appreciated that in addition to tumor cell-intrinsic survival and growth 54 programs, the microenvironment is crucial in supporting metastases formation 1-3 .

123
Initial data analysis indicated that fibroblasts isolated from lungs with macrometastases 124 (macrometastasis-associated fibroblasts-MAF) were strikingly different from NLF as well as 125 from fibroblasts isolated from lungs with micrometastases (micrometastasis-associated 126 fibroblasts-MIF) ( Figure 1H,I, Figure 1 -figure supplement 1). Notably, since fibroblasts were 127 isolated from entire lungs, rather than from specific metastatic lesions, the MIF fraction 128 contained a mixture of normal, non-metastasis-associated fibroblasts as well as metastasis-129 associated fibroblasts. As a result, initial data analysis did not reveal significant differences 130 between NLF and MIF. Thus, metastasis-associated fibroblasts are not only functionally 131 activated but also transcriptionally reprogrammed.

133
Transcriptome profiling of metastasis-associated fibroblasts reveals dynamic stage-134 specific changes in gene expression.

135
In light of these initial results, we next analyzed the genes that are differentially expressed 136 between MAF and NLF. We selected upregulated and downregulated genes based on fold 137 change of |2|. Expectedly, hierarchical clustering based on these genes revealed that the 138 MAF group clustered separately from NLF and MIF (Figure 2A). To better characterize the 139 trajectory of changes in fibroblasts during metastases formation, we next compared the 140 expression of genes that were differentially expressed between MAF and NLF to their 141 expression in the MIF population. Interestingly, we found that the expression pattern in MIF 142 was distinct from both the MAF and the NLF gene expression, including genes that had 143 opposite changes in MAF vs. MIF, suggesting that they are activating a distinct 144 transcriptional program ( Figure 2B).

145
We therefore analyzed the differentially expressed genes in the MIF fraction separately.

146
Since the detectible changes in micrometastases were more subtle than the changes 147 detected in the macrometastases group, we selected these genes based on a fold change of 148 |1.5|, to better differentiate the MIF group from NLF. Indeed, hierarchical clustering based on 149 these differentially expressed genes confirmed that the MIF group clustered separately from 150 both NLF and MAF ( Figure 2C). Next, we selected a group of genes based on their 151 differential expression between the MAF and MIF groups (FC>|2|). The combination of these 152 yielded a total of 897 genes that were differentially expressed in MIF vs. NLF, MAF vs. NLF

182
Interestingly, we found that gene expression signatures in fibroblasts isolated from the micro-183 metastatic stage were highly and specifically enriched for functions related to cellular 184 response to stress, including Hsf1 activation, heat shock response and response to unfolded 185 protein (Supplementary File 1). Upregulated genes in MIF that were related to stress and 186 protein folding included several heat shock proteins: Hspa8, Hsp90aa1, Hspd1, Hspe1 and 187 others ( Figure 3C). Of note, detailed analysis of specific gene expression showed that while 188 the stress response pathway was not significantly enriched in MAF, genes from the stress 189 response pathway were elevated in MAF compared to normal fibroblasts, but not compared 190 to MIF ( Figure 3C). ECM remodeling terms were enriched in both MIF and MAF ( Figure 3B),

191
indicating the central importance of ECM modifications in facilitating metastasis. Notably, 192 while ECM remodeling was operative throughout the metastatic process, the specific genes 193 related to ECM remodeling in the different metastatic stages were distinct ( Figure 3D).

194
Gene expression signatures in fibroblasts isolated from macrometastases were highly 195 enriched for inflammation-related pathways ( Figure 3B, Supplementary File 1). Indeed, 196 analysis of enriched pathways revealed that genes related to inflammation including many 197 chemokines and cytokines were upregulated specifically in MAF ( Figure 3E). To validate 198 these findings, we isolated fibroblasts from additional cohorts of mice. We performed qRT-199 PCR to test the expression of key genes from identified pathways (stress response, ECM

210
Taken together, these findings imply that metastasis-associated fibroblasts assume distinct 211 functional roles during the process of lung metastasis.

212
Encouraged by these findings, we next set out to obtain further insights on functional 213 pathways that were modified in fibroblasts isolated from different metastatic stages. To that 214 end, we performed Gene Set Enrichment Analysis (GSEA) 39 . We focused our analysis on 215 the H collection: Hallmark gene sets that summarize specific well-defined biological states or 216 processes based on multiple datasets 40 . Similar to the results obtained in our previous 217 analyses, we found that functions related to inflammatory responses, including TNF and IL-

231
To further characterize the regulatory nodes that govern the transcriptional changes in 232 fibroblasts, we hypothesized that these changes may be driven by transcription factors (TFs) 233 related to the pathways that were identified by the pathway and GSEA analyses (Figure 3).

234
Analysis of TFs terms within the results identified five candidate transcription factors (TFs) that were enriched in at least one analysis and in at least one metastatic stage: Hif1a, Hsf1,

237
We next examined the number of different comparisons in which each TF was enriched. We

242
To rank these TFs, we performed knowledge-based multiple analyses examining their  analyzed the correlation of the metastasis-associated gene network with each candidate TF 265 using the VarElect tool 43 . This tool enables prioritization of genes related to a specific query 266 term by using a direct and indirect relatedness score. We analyzed the scores of the stage-267 specific signature genes with each candidate TF, and the number of directly related genes.

268
The TFs were ranked based on the number and average score for the directly related genes, 269 and the average score of the indirectly related genes. In agreement with previous analyses,

270
Myc had the highest number of connections and the highest average score for both directly 271 and indirectly related genes in all comparisons ( Figure 4A, pink, Figure 4C). To consolidate 272 these comprehensive gene network analyses, we performed a comparative analysis on the 273 TF bioinformatics measurements listed in Figure 4A. The results indicated that Myc achieved 274 significantly higher scores than all other TFs in all three gene signatures ( Figure 4D).

275
Since the changes in transcriptome were associated with multiple TFs, we further asked

292
To validate the ranking results, we analyzed by qRT-PCR the expression of Myc in 293 fibroblasts isolated from normal lungs, or from lungs with micro-and macrometastases.

294
Analysis of the results indicated that Myc is significantly upregulated in macrometastases-295 associated fibroblasts ( Figure 5A). In addition, we assessed the expression of central Myc 296 targets that we found to be upregulated in metastasis-associated fibroblasts, including 297 Hspe1, Hsp90aa1, Odc1 and Fosl1 46,47 . The results indicated that these Myc targets were 298 upregulated in fibroblasts isolated from lungs with metastases ( Figure 5B). qRT-PCR results

299
of Myc target genes further confirmed that the stress response-related genes Hsp90aa1 and

408
In addition to stress response, our findings indicated that ECM remodeling is a central task of

459
The authors would like to thank Dr. Ran Elkon for his help with data analysis.

462
The authors declare no conflict of interests.

466
All experiments were performed using 6-8 weeks old female mice, unless otherwise stated.

467
All experiments involving animals were approved by the Tel Aviv University Institutional

540
Sections were fixed with 4% PFA for 5 min, permeabilized by 0.2% Triton for 20 min and 541 fixed with NBF as described above. Antigen retrieval was performed using citrate buffer (pH 542 6.0). Slides were blocked with 1% BSA, 5% normal goat serum in 0.2% PBST for 1h and

625
Results were considered significant with a p-value<0.01, q-value<0.05 and a coverage >3%.

626
To increase the specificity of the enriched terms, we compared the relative overlap and the

700
All experiments represent at least 3 separate biological repeats, unless otherwise stated.

806
Data are represented as mean ± SD, n=3.