Transcriptomic profiles conducive to immune-mediated tumor rejection in human breast cancer skin metastases treated with Imiquimod

Imiquimod is a topical toll-like-receptor-7 agonist currently used for treating basal cell carcinoma. Recently, imiquimod has demonstrated tumor regression in melanoma and breast cancer skin metastases. However, the molecular perturbations induced by imiquimod in breast cancer metastases have not been previously characterized. Here, we describe transcriptomic profiles associated with responsiveness to imiquimod in breast cancer skin metastases. Baseline and post-treatment tumor samples from patients treated with imiquimod in a clinical trial were profiled using Nanostring technology. Through an integrative analytic pipeline, we showed that tumors from patients who achieved a durable clinical response displayed a permissive microenvironment, substantiated by the upregulation of transcripts encoding for molecules involved in leukocyte adhesion and migration, cytotoxic functions, and antigen presentation. In responding patients, Imiquimod triggered a strong T-helper-1 (Th-1)/cytotoxic immune response, characterized by the coordinated upregulation of Th-1 chemokines, migration of Th-1 and cytotoxic T cells into the tumor, and activation of immune-effector functions, ultimately mediating tumor destruction. In conclusion, we have shown that topical imiquimod can induce a robust immune response in breast cancer metastases, and this response is more likely to occur in tumors with a pre-activated microenvironment. In this setting, imiquimod could be utilized in combination with other targeted immunotherapies to increase therapeutic efficacy.

In melanoma, several case reports have shown regression of cutaneous metastases with imiquimod either alone 18,19 or in combination with other therapies [20][21][22] . This tumor regression was associated with an increase in T cell infiltrate and upregulation of ICR genes such as CCL5, CXCL9, CXCL10 important for homing of cytotoxic T cells as well as markers of dendritic cell (CD80, CD86) and T cell activation (CD69) 23 .
Activation of the adaptive immune response has been shown to be associated with better prognosis in breast cancer. An increase in the number of tumor infiltrating lymphocytes is associated with greater likelihood of complete pathological response after chemotherapy across breast cancer subtypes, and reduced risk of relapse and death particularly in triple negative breast cancers [24][25][26][27][28][29][30][31] . Interestingly, most of the prognostic and predictive transcriptomic classifiers described so far in breast cancer include ICR genes. Type I immunity is believed to be especially important for immunotherapy, with CD8+ T cells inducing cytotoxicity and CD4+ T helper cells inducing cytokine secretion and promotion of antigen presentation of tumor 32 . In a recent trial, topical imiquimod in combination with intravenous nab-paclitaxel resulted in an antitumor response of 72% in the treatment of mostly hormone receptor negative breast cancer cutaneous metastases 33 suggesting that activation of the immune system may play an important role in the treatment of breast cancer skin metastases. Imiquimod may play an even more important role in hormone receptor positive breast cancer, since it is the subtype with the least amount of CD8+ infiltrating T cells 34 .
We have previously shown that topical imiquimod, as a single agent, can induce an antitumor response in 20% of skin metastases of mainly hormone receptor positive breast cancers and can promote an immunogenic tumor microenvironment 35 . An increase in CD4+ and CD8+ infiltrating T cells after treatment was seen in one of the responders but not in the other, and cytokine measurement in tumor supernatant was variable, suggesting the need for a more sensitive assay to characterize the effects of imiquimod on the tumor microenvironment in breast cancer skin metastases. Furthermore, in two additional patients who achieved long-term benefit, an antitumor antigen response was induced by imiquimod (in-situ vaccination effect) and subsequently expanded by endocrine therapy leading to durable complete remissions. 36 Here, by using an integrative analysis we describe transcriptomic profiles associated with responsiveness to imiquimod treatment. This is the first study to characterize the transcriptomic changes induced by imiquimod in breast cancer skin metastases.

Methods
Patients. Ten patients were enrolled and treated with imiquimod for eight weeks as previously described 35 .
The clinical trial was approved by the New York University Institutional Review Board. All research was performed in accordance with the New York University Institutional Review Board guidelines and regulations, a written informed consent was obtained from all patients. Same-site tumor biopsies were obtained at baseline and after 8 weeks of imiquimod treatment. Paraffin embedded tumor tissue was available from 8/10 patients for this study, as samples from two patients had insufficient quantity for further analysis. Two of the patients had stable disease during the initial study and were found to have a systemic complete clinical response after subsequent treatment with fulvestrant after study completion (including the treated skin metastases). On follow up, these two patients also had disease remission for two years. We initially sought to characterize the tumor microenvironment of these two patients due to their unusual complete response on hormone therapy and long duration of disease remission and labeled these patients as complete responders (CR). An additional patient had a local partial tumor response (greater than 50% tumor shrinkage) after eight weeks of imiquimod treatment and was labeled as a partial responder (PR). Patients with CR and PR were considered as having derived clinical benefit (CB). Five of the eight patients did not have a clinical response and were defined as non-responders (NR). Gene Expression. Paraffin blocks were carefully evaluated for tumor content and cut into seven sections of 10 um thickness for RNA extraction. Samples were deparaffinized with xylene and extracted using the RNeasy FFPE kit (Qiagen #73504), according to manufacturer's instructions. Extracted FFPE RNA quality and quantity were analyzed on an agilent Bioanalyzer 2100 using a nano chip. Smear analysis was performed to determine the percent of RNA greater than 300nt for optimized hybridization according to Nanostring's protocol. Adjusted inputs for each RNA sample were calculated to input 100 ng of RNA greater than 300nt. RNA was hybridized using the Nanostring nCounter ® Human v1.1 PanCancer Immune Profiling Panel (770 transcripts) according to the manufacturer's protocol. Hybridizations were processed on the nCounter Prep Station, and prepped cartridges were scanned on the Nanostring Digital Analyzer using 280 field of view counts.
Immunohistochemistry. The evaluation of tumor-infiltrating lymphocytes (TIL) density was performed by Immunohistochemistry (IHC) on paraffin embedded tumor tissues as previously described 35 and correlated with clinical response. Tumor sections (thickness 4 μm) were deparaffinized and rinsed in distilled water. Heat induced epitope retrieval was carried out in 10 mmol/L citrate buffer. CD3, CD4, CD8 (Ventana Medical Systems), and Forkhead Box Protein P3 antibody (FoxP3, Ebiosciences) antibodies were used. Detection was carried out on a NEXes instrument (Ventana Medical Systems) using the manufacturer's reagent buffer and detection kits. After washing in distilled water, slides were counterstained with hematoxylin, dehydrated, and mounted with permanent media. Appropriate positive and negative controls were included with the study sections. IHC-positive cells were counted manually in 5 representative high-power fields (HPF, 400), to derive the average number per HPF, by a pathologist blinded to the response. Statistical Analysis. Data were normalized using the housekeeping genes present in the panel and by applying negative control subtraction (nSover 2.6 package). Quantile normalization using preprocessCore v1. 38 and Log 2 transformation was applied to the expression matrix. PCA plots were generated using scatterplot3d v0.3. Consensus clustering based on the ICR genes was performed using consensusclusterPlus v1.40 (maxK = 7, 5000 repetitions, Ward D2 as interlinking). Based on such clustering, samples were classified as ICR High, ICR Medium, and ICR Low, with ICR High and ICR Low samples as the highest and lowest expression of the ICR z-score, respectively.
Singe sample gene set enrichment (ssGSEA) 37,38 using GSVA package v1.24.2 was applied to calculate enrichment scores for leukocyte populations using cell-specific transcripts included in the Nanostring panel (Supplementary Table S1) as well as the enrichment score for the ICR genes. The fold changes of these scores were expressed as the ratio of Anti Log e of the mean enrichment scores and were plotted with their confidence interval using forestplot v1.7.2. Differentially expressed genes (DEGs) or differences in the enrichment scores were calculated by paired or unpaired t-test for the post-vs pre-treatment and CR vs NR analysis respectively using 2-tailed p value of 0.05 as significance cut off. Fold change was expressed by the Anti-Log 2 of the Log 2 fold change, which is the difference between the mean of the Log 2 values of the two categories in comparison. Enrichment analysis of the DEGs was performed using Nanostring annotations (immune response categories, Supplementary  Table S2) specific for this panel. In all comparisons of response status, only the CR and NR samples were included as these samples represent tumors with two opposite biologic behaviors. PR data were however displayed in the heatmaps and boxplots in comparison with CR and NR samples. Post-vs pre-treatment comparisons were performed by including all samples. Boxplots for single gene expression and ICR score were plotted using ggplot2 v 2.2.1. Dotted heatmaps were used to represent the ssGSEA enrichment z-scores and the percentage of genes with an absolute fold change of more than 1.5 for each immune cell signature. These plots were generated using the ComplexHeatmap v1.14 package. Immunohistochemistry data were compared by Wilcoxon signed-rank test.

Results
Patient characteristics. Paired tissues from eight patients with breast cancer skin metastases enrolled on a clinical trial of imiquimod therapy were included in this study. Patient demographics are summarized in Table 1. The majority of the women were post-menopausal (62.5%), all had invasive ductal carcinoma, and the majority of tumors were high grade (87.5%) and hormone receptor positive (87.5%), 50% were HER2 positive. Patients were designated as CR (n = 2), PR (n = 1), and NR (n = 5) as described above.

Transcriptome changes induced by imiquimod treatment and baseline predictive markers.
Principal component analysis (PCA) based on all the transcripts is shown in Fig. 1A. Using a PC1 arbitrary cut-off of −10 3/3 tumors from patients responding to the treatment but only 1/5 tumors from non-responding patients were found to cluster. The number of overlapping genes among the comparisons used for the detection of the DEGs (e.g., post vs pre-treatment, pre-treatment CR vs NR, and post-treatment CR vs NR) are shown in Fig. 1B. Comparing all patients post-versus pre-treatment, 53 DEGs were identified ( Table 2, Fig. 1B,C). Several genes related to innate immunity were upregulated in the post-treatment group. These included genes related to lymphocyte activation (SH2D1A, SLAMF1), lymphocyte adhesion (SELL), cytolysis (PRF1), NK cell activation (KLRC1), T cell (CD3G, CD347) and cytotoxic T cell (CD8A) function, and antigen presentation (HLA-C). Although some genes were downregulated, overall, pathway enrichment analysis revealed a positive modulation of the immune response underlined by the activation of cytotoxic mechanisms ( Fig. 1D  www.nature.com/scientificreports www.nature.com/scientificreports/   Table 2. DEG post vs pre-treatment. 53 Significant (p < 0.05) DEGs comparing nanostring immune gene expression between post-and pre-treatment samples using paired t-test including non responder (NR), complete responder (CR) and partial responder (PR) samples. Ordered from high to low fold change (Fc = 2^(mean log2 post treatment expression − mean log2 pre treatment expression)).
In baseline samples, prior to treatment with imiquimod, 45 genes were differentially expressed in CR vs NR samples (Table 3, Fig. 1B). Interestingly, the PR clustered between CR and NR (Fig. 1C). Upregulated genes included leukocyte adhesion and migration molecules (ITGAX, ITGAM, ITGB2, and FUY7), chemokine  www.nature.com/scientificreports www.nature.com/scientificreports/ receptor (CX3CR1), T-cell cytotoxic markers (CD8A, and GZMK), and genes associated with antigen presentation (HLA-DMA, HLA-DBP1, and HLA-DPA1). Enrichment analysis corroborated the existence of a permissive microenvironment in lesions from patients with clinical benefit from imiquimod, substantiated by the activation of the antigen processing pathway, which was accompanied by the downregulation of some cytokines such as IFN lambda 1 (IFNL1) IL34, and IL23, which is involved in Th-17 signaling (Fig. 1D).
Applying multiple hypothesis correction of a 0.25 false discovery rate (q values provided in Tables 2, 3

Responders to imiquimod treatment had higher levels of T cells and cytotoxic cells. Single sam-
ple gene set enrichment analysis (ssGSEA) was used to estimate the variation in leukocyte subpopulations in each sample. Post-vs pre-treatment comparison revealed that overall, imiquimod consistently enhanced T cell (p = 0.011) and cytotoxic T cell (p = 0.017) infiltration ( Fig. 2A). Before treatment, CR patients had a lower level of T central-memory (Tcm) cells (p = 0.038) and a higher level of macrophages (p = 0.0057), T cell (p = 0.039) and NKCD56 bright cells (p = 0.012) as compared to NR patients. The relative depletion of Tcm cells could be the result of a proportional enrichment of activated T-cell after treatment. Thus, CR patients displayed a baseline pre-activated immune microenvironment (Fig. 2B).
In CR patients, there was a strong activation of the cytotoxic response (p = 0.012) accompanied by upregulation of the ICR genes ( Fig. 2C) (p = 0.023) and by a switch from immature dendritic cells (iDC) (p = 0.0228) to activated dendritic cells (aDCs) (p = 0.096) in post-treatment compared to baseline biopsy. At the individual sample level, the magnitude of the cytotoxic response was maximal in responding patients, with similar levels among the two CR and the PR (Fig. 3A,B).
A nonsignificant trend for higher CD3, CD8, CD4, and FOXP3 positive cells, as evaluated by immunohistochemistry, was observed in CR vs NR pre-treatment samples and was also present after treatment. However, among the responding patients, IHC evaluation detected increase of TILs (consistent across all the markers) only in the PR patient (Supplementary Figure 1), who had a local response at the time of the biopsy. This observation suggests that molecular profiling might be more sensitive than IHC in detecting changes in the functional orientations of the tumor microenvironment preceding the clinical response.
Responders to imiquimod treatment displayed consistent upregulation of the ICR score. We analyzed the immunologic constant of rejection (ICR) score (ie., the mean log2 expression of the ICR genes, Fig. 4A,B) across samples, and compared its expression in CR versus NR, pre-and post-treatment. The ICR score was significantly higher in the post-treatment CR vs NR samples (p = 0.021) while it did not differ among the two categories at baseline (Fig. 4C). Interestingly, the PR sample clustered with the CR samples. The ICR genes with significantly higher levels in CR vs NR samples after treatment are shown in Fig. 4D. Most of the ICR genes had similar and consistent trends (Supplementary Figure 2). An exception was represented by FOXP3, which displayed a non significant lower level in post-treatment CR vs NR samples.

Discussion
This is the first study to characterize the transcriptional changes induced by imiquimod in breast cancer skin metastases. The CR samples showed a pre-activated immune microenvironment at baseline, substantiated by the enrichment of genes related to antigen processing and a proportional enrichment of T cells, NKCD56 bright cells and macrophages. Macrophages have often been shown to be a negative prognostic factor, although M1 macrophages, which sustain anti-tumor immune response, have been associated with favorable prognosis 39,40 .  Table 4. DEG CR vs NR (post-treatment). 57 Significant (p < 0.05) DEGs comparing nanostring immune gene expression in post treatment samples between CR and NR using unpaired t-test (PR sample was excluded). Ordered from high to low fold change (Fc = 2^(mean log2 CR − mean log2 NR expression)).  www.nature.com/scientificreports www.nature.com/scientificreports/ Only a limited proportion of DEG overlaps between pre-treatment CR vs NR, post-treatment CR vs NR, and post-vs pre-treatment lesions. The lack of overlap among DEG post-vs pre-treatment with the ones detected in the post-treatment CR vs NR comparison might be the result of the dilution of the signal in the post-vs pre-treatment comparison due to the inclusion of the NR samples. The differences observed in pre-and post-treatment CR vs NR samples might be the consequence of a differential functional activation of the immune response at these time points. Baseline CR samples are dominated by a sub-inflamed microenvironment typified by the upregulation of antigen-processing related genes, while post-treatment differences are characterized by the www.nature.com/scientificreports www.nature.com/scientificreports/ coordinated activation of cytotoxic mechanisms conducive to immune-mediated rejection. Nevertheless, 40% (3/7) of the upregulated genes overlapping between at least two comparisons reflect the activation of a cytotoxic response (CD8A, GZMK, KLRK1) while the remaining ones are related to immune modulation. This observation suggests that the full activation of a pre-existing, yet incomplete, cytotoxic response is critical to induce immune-mediated tumor rejection in this setting.
LAG-3 was identified as one of the DEGs upregulated with a 6.6 fold change in the CR group post-treatment. LAG-3 is a receptor expressed by activated T lymphocytes and is considered a marker of T cell exhaustion leading to reduced effector function and expression of inhibitory receptors 43 . Blockade of LAG-3 is hypothesized to synergize with other immune checkpoint inhibitors such as anti-PD-1 44 . Phase II data of imiquimod with nab-paclitaxel also showed that patients with a poor clinical response had higher circulating levels of PD-1+ T cells suggesting that combination therapy with anti-PD-1 and/or anti-LAG-3 may increase treatment responses 33 . Previous clinical trial data, of mostly estrogen receptor positive metastatic breast cancer, suggest that when combined with chemotherapy, a LAG-3 inhibitor in metastatic breast cancer can increase the number of activated monocytes, dendritic cells, NK cells and CD8 T cells with promising response rates 45 . Our data suggests that LAG-3 may play an important role in breast cancer skin metastases as well and LAG-3 inhibition and PD-1 inhibition should be considered in combination with imiquimod as possible treatment modalities.
Interestingly, the PR sample had the same immune gene expression signature post-treatment (clustering together with the 2 CR samples) as the complete responders (Fig. 1C). Our subsequent long-term analyses of tumor antigen-specific immune responses point to the important role of imiquimod treatment inducing the adaptive response, which then got augmented during endocrine therapy in both patients as they continued to derive clinical benefit and entered a long-lasting CR 36 . This likely explains the unexpected and durable complete response from fulvestrant since a previous large phase 3 study of fulvestrant monotherapy demonstrated that a complete response is very rare 46 . This further suggests that the genes upregulated by imiquimod are associated with long term immunogenic effects which can improve the efficacy of subsequent therapies. Whereas upregulation of ICR genes has been previously correlated with better response to chemotherapy and trastuzumab 29,47 , our data suggests that activation of these pathways may also improve response to subsequent hormonal therapy.
Limitations of our study are the small number of patients and the short duration of treatment. The delayed CR in two patients who displayed immune activation within the metastases after eight weeks of treatment with imiquimod suggests the possibility that a longer treatment might increase response rates. Genes differentially expressed detected by univariate analysis are subject to false discovery. The possibility to incur into a Type I error (i.e., the rejection of a true null hypothesis) is partially attenuated by the use of the integrative pipeline presented here. By combining deconvolution and enrichment analyses, more permissive cut-off p values can be used to detect biological relevant phenomena. Nevertheless, this dataset should be considered as a discovery set, and further validation using orthogonal platforms as deep phenotyping and proteomic approaches are warranted.
The strength of this study is the use of gene array analysis in paired samples derived from a prospective clinical trial with long-term follow-up of patients. This sensitive technique for characterizing immune changes induced by imiquimod treatment was able to detect at 8 weeks intratumoral changes associated with not only immediate but also future clinical benefit, which could not be derived from histologic analysis 35 .
Metastatic breast cancer, especially with cutaneous metastases, remains very difficult to treat. Long-term clinical response to therapy in metastatic cancer is rare. Promoting immunogenicity in the tumor microenvironment may enhance efficacy of other therapies such as chemotherapy, endocrine therapy, targeted therapies or even newer immunotherapeutics. Imiquimod has been shown to be a potent activator of the innate immune system, which then supports adaptive immune response in cancer, and we demonstrate that this is the case in metastatic breast cancer as well. Encouragingly, although hormone receptor positive breast cancer appears to be the least immunogenic type of breast cancer 34 , we have shown that imiquimod can induce a robust anti-tumor immune response leading to clinical response in the three patients with hormone receptor positive disease described here. This is in line with a recent report of successful immunotherapy in a patient with metastatic hormone receptor positive breast cancer 48 . Future trials are needed to investigate the use of imiquimod in combination with other therapies.

Data Availability
The normalized Nanostring data used in the current study is available as Supplementary File 1.