Research ArticleCANCER

Tumor immunological phenotype signature-based high-throughput screening for the discovery of combination immunotherapy compounds

See allHide authors and affiliations

Science Advances  22 Jan 2021:
Vol. 7, no. 4, eabd7851
DOI: 10.1126/sciadv.abd7851

Abstract

Combination immunotherapy is promising to overcome the limited objective response rates of immune checkpoint blockade (ICB) therapy. Here, a tumor immunological phenotype (TIP) gene signature and high-throughput sequencing–based high-throughput screening (HTS2) were combined to identify combination immunotherapy compounds. We firstly defined a TIP gene signature distinguishing “cold” tumors from “hot” tumors. After screening thousands of compounds, we identified that aurora kinase inhibitors (AKIs) could reprogram the expression pattern of TIP genes in triple-negative breast cancer (TNBC) cells. AKIs treatments up-regulate expression of chemokine genes CXCL10 and CXCL11 through inhibiting aurora kinase A (AURKA)–signal transducer and activator of transcription 3 (STAT3) signaling pathway, which promotes effective T cells infiltrating into tumor microenvironment and improves anti-programmed cell death 1 (PD-1) efficacy in preclinical models. Our study established a novel strategy to discover combination immunotherapy compounds and suggested the therapeutic potential of combining AKIs with ICB for the treatment of TNBC.

INTRODUCTION

Although ICB therapy is now understood as a powerful therapeutic approach for treating and even curing patients with cancer, there is a well-recognized problem with low objective response rates (ORRs) for a majority of solid tumors (1, 2). Unfortunately, many checkpoint blockade immunotherapies simply do not work for large proportions of patients, with ORRs being as low as 10% for common highly aggressive cancers such as TNBC. There is an emerging consensus that these low ORRs are due, at least partially, to poor infiltration of effective T cells into the tumor microenvironment (3). As low ORRs have brought tumor heterogeneity into clearer focus in recent years, researchers have been paying increasing attention to the so-called TIP (4); for example, tumors are now frequently categorized as either “cold” (noninflamed) or “hot” (inflamed) on the basis of the extent of immune cell infiltration.

Genetically, it is known that immune cell infiltration is regulated by tumor-intrinsic oncogenic and epigenetic pathways. Cell biology studies have also shown that immune suppressive cells negatively affect T cell infiltration. Extensive clinical work has identified particular genes that are strongly associated with either cold tumors or hot tumors; very briefly, widely recognized hot tumor genes include, among many others, T helper 1 (TH1)–type chemokines such as CXCL9, CXCL10, and CXCL11 that promote the infiltration of effector immune cells (e.g., CD8+ T cells, TH1 cells, TH17 cells, and antigen-presenting cells). In this context, one can understand the desirability of identifying small molecules that could convert immunologically cold tumors into hot tumors and might markedly increase ORRs for immunotherapies (5). Foundational studies demonstrated that a small molecule can epigenetically alter the expression of genes known to affect immune cell infiltration (TH1 chemokines) and thereby confers a therapeutic effect against ovarian and colon cancer cells (6, 7).

Combination therapy is the key to overcome the low ORRs of checkpoint blockade immunotherapy. One purpose of combination strategy is to improve therapeutic efficacy by overcoming the immunosuppressive microenvironment of tumors. Known combinations of the checkpoint inhibitors explored to date include those with (i) other immune checkpoint inhibitors (8), (ii) chemotherapy (9), and (iii) targeted therapies (10, 11). However, almost all these combination agents are selected on the basis of known targets or known drugs. It is hard to discover unknown targets or compounds for the immunotherapy combinations, partially because of the absence of effective and high-throughput screening approaches.

The gene signatures of tissues or cells have been widely used to indicate cancer metastasis (12) and subtype cancers (13), to estimate clinical outcome of patients (14) and to evaluate tumor purity, and to infer the abundance of immune cell types (15), as well as to predict cancer immunotherapy response in recent research (16). In addition, gene expression profiling of cellular perturbations has also been extensively used to discover anticancer drugs (17). Here, specifically, we generated a gene signature, which could distinguish the immunosuppressive microenvironment of cold tumor to that of hot tumor, and used the expression changes of this gene signature as the read-out for high throughput of drug screenings. We decided to use high-throughput sequencing–based high-throughput screening (HTS2) (18) for the discovery of immunotherapy combination agents. HTS2 takes advantage of the next-generation sequencing technologies and has been successfully applied for the discovery of novel compounds against prostate cancer (18) and against breast cancer lung metastasis (19).

Given its extremely low ORRs, we were particularly interested in TNBC, an aggressive breast cancer subtype defined by a lack of expression of estrogen receptor, progesterone receptor, and human epidermal growth factor receptor 2 (3, 20). About 15% of all newly identified breast cancers are diagnosed as TNBC, which is associated with an increased risk of distant metastasis and death compared with other breast cancer subtypes (21). This poor prognosis is further compounded by the fact that patients with TNBC are insensitive to most available therapies, including immunotherapies. Clearly, effective therapies for TNBC are urgently needed. The aim of the present study was to explore the possibility of using TIP gene signature and HTS2 as a new approach to discover potential immunotherapy combination agents for the treatment of TNBC.

RESULTS

TIP gene signature is associated with the immune state of tumor

Our HTS2 strategy requires a defined gene expression signature that reflects the functional process tested by the screen (here, immune cell infiltration). To initially identify candidate genes for our targeted immunologic phenotype, we queried >30,000 relevant literatures using information extraction (IE), a text-mining approach (22). The top 15 genes identified in the IE analysis were initially selected as our “TIP genes” (Fig. 1A). This set included 12 hot tumor–related genes (CXCL9, CXCL10, CXCL11, CXCR3, CD3, CD4, CD8a, CD8b, CD274, PDCD1, CXCR4, and CCL5) and 3 cold tumor–related genes (CXCL1, CXCL2, and CCL20). Encouragingly, by the analysis of a published melanoma dataset including 266 metastatic human cutaneous melanoma samples (23), we found that these hot tumor– or cold tumor–related genes we selected are highly expressed in reported hot or cold melanoma tumors, respectively (Fig. 1B). Moreover, the expression of our selected TIP signature genes can be used to distinguish hot with cold tumors with 89.7% accuracy (Fig. 1C) and good performance (Fig. 1D).

Fig. 1 The expression pattern of TIP signature genes is associated with TIPs.

(A) TIP signature that contains 15 immune genes. Red column, highly expressed genes. NK, natural killer. (B) Expression pattern of TIP signature genes within the hot and cold melanoma cohorts. (C) The TIP signature could distinguish hot and cold tumors with 89.7% accuracy in a melanoma dataset. (D) The receiver operating characteristic curve analysis was applied to evaluate the accuracy of the signature genes for tumor classification in a melanoma dataset (23). AUC, area under curve. (E) The correlation analysis of the expression level of six TIP signature genes with CD8+ T cell infiltration level in breast invasive carcinoma (METABRIC n = 1980); RSEM, RNA sequencing (RNA-seq) by expectation-maximization.

We next checked whether the expression of these TIP signature genes was actually related to the immune state of tumors by examining correlations between their expression levels and the expression levels of known immune cells marker genes and the infiltration level of CD8+ T cell in a breast tumor dataset [METABRIC (24)]. Known markers for the presence of effective T cell in tumors include CD8, CD4, and CD3e (25). By analyzing 1980 invasive breast tumors from a dataset [METABRIC (24)], we found that the expression of at least six of our TIP signature genes was notably associated with these three marker genes (fig. S1, A and B) and the CD8+ T cell infiltration level (Fig. 1E).

TIP gene signature is associated with cancer immunotherapy responses

Next, we wondered how the TIP gene signature performed on predicting ICB response. To evaluate this, we collected four ICB therapy datasets including three with patients with melanoma (2628) and one with patients with glioblastoma (29) and checked the correlation of the expression level of TIP gene signature and the response of ICB therapy patients. As expected, in these four ICB clinical cohorts, hot tumor–related genes were consistently higher in responders than nonresponders and vice versa (Fig. 2, A and B). In addition, patients in hot tumor–related genes highly expressed group have better prognosis (Fig. 2C).

Fig. 2 The expression pattern of TIP signature genes predicts the response of immunotherapy.

(A) Averaged expression level of TIP-associated gene in pretreatment tumors of ICB responders and nonresponders over four cohorts (2629), represented separately (left) and merged (right). Within each group, the scattered dots represent the averaged z score of hot tumor–related genes of each sample, and the thick line represents the median value; the bottom and top of the boxes are the 25th and 75th percentiles (interquartile range). Whiskers encompass 1.5 times the interquartile range. (B) Heatmap showing the expression of TIP-associated gene that discriminates between anti–PD-1 responders and nonresponders. (C) Kaplan-Meier plots of overall survival (OS) for 42 patients with melanoma who accepted ICB therapy (28) and were grouped with hot tumor–related genes highly or lowly expressed. (D and E) AUC for signatures in predicting ICB response in seven datasets (27, 28, 3135), represented merged (D) and separated (E). AUC in (D) is averaged AUC among seven datasets; the performance of a random predictor (AUC, 0.5) is represented by the dashed line. NA, the data are not available. TMB, tumor mutation burden; MSI, microsatellite instability signature.

We further compared TIP gene signature with eight widely used ICB response biomarkers including IFNG, CD8, programmed cell death 1 ligand 1 (PD-L1), tumor mutation burden, microsatellite instability signature, tumor immune dysfunction and exclusion (TIDE), T cell clonality, and B cell clonality with biomarker evaluation tool of TIDE (30). Among all candidate biomarkers, we found that TIP gene signature performed as the best predictor over seven datasets with both anti–PD-1 and anti-CTLA4 therapies (Fig. 2, D and E) (27, 28, 3135).

TIP gene signature is associated with clinical outcomes of patients with cancer

We were also interested in whether our selected TIP signature genes were associated with the clinical outcomes of patients with cancer. By calculating the log-rank P value (Fig. 3A) and hazard ratio (Fig. 3B) of TIP signature genes on overall survival (OS) over 28 The Cancer Genome Atlas (TCGA) cancer types with multivariable Cox proportional hazard model, we found that the expression of TIP signature genes were correlated with the OS in most of the cancer types (Fig. 3, A and B) and with breast cancer ranked at the first (Fig. 3B). We then zoomed our analysis in patients with breast cancer.

Fig. 3 The TIP gene signature is correlated with the clinical outcome of patients with cancer.

(A and B) The impact of TIP signature genes on OS in patients over 28 cancer types, dotplot with log-rank P value (A) and forest plot of Cox proportional hazard ratio and 95% confidential interval (CI) of OS (B) were showed. (C and D) The genomic alterations (C) and the percentage of patients with each genetic alteration type (D) for TIP signature genes in the TCGA breast cancer dataset (58). (E) Fisher’s exact test of genetic deletion (or mRNA down-regulation or missense mutation) and amplification (or up-regulation) for hot tumor– and cold tumor–related genes. (F and G) Impact of TIP signature genes on OS in patients with breast cancer; log-rank test was used. (H) Top: Ranked risk score of 1980 patients with breast cancer evaluated by the expression and the risk coefficient of TIP signature genes. Middle: OS distribution of the 1980 patients ranked according to the risk score (from the top). Bottom: Heatmap of the expression pattern of each of the TIP signature gene in the 1980 patients.

The analysis of the breast cancer dataset [METABRIC (24)] revealed that the cold tumor genes of this signature were generally amplified and/or up-regulated, while the hot tumor genes were very frequently deleted and/or down-regulated (Fig. 3, C to E). Kaplan-Meier analysis indicated that the individual expression of the cold tumor genes CXCL1 and CXCL2 in this signature were in a trend of negatively correlated with OS (Fig. 3F), whereas the expression of the hot tumor genes CXCL10, CXCL11, CCL5, and CXCR3 were positively correlated with OS (Fig. 3G).

We also performed risk analysis for individual patients by using the expression level and coefficient of individual TIP signature genes in breast tumor dataset [METABRIC (24)]. These patients were classified into either low- or high-risk group (Fig. 3H), and our analysis indicated that cold tumor genes of this signature such as CXCL2, CXCL1, and CCL20 were up-regulated in the high-risk group and that hot tumor genes such as CXCL9, CXCL10, and CXCL11 were up-regulated in the low-risk group. Collectively, these results imply that our selected TIP signature genes are meaningfully associated both with the effective T cell amount in tumors, with the ICB response and with the prognosis of patients with breast cancer.

HTS2 screening identified aurora kinase inhibitors that reprogram the expression of TIP signature genes in TNBC cells

Given the significant correlation of TIP gene signature with the effective T cell amount in tumors, with the ICB response, and with the prognosis of patients with breast cancer, we decided to perform a high-throughput drug screening to identify some compounds to revert the expression pattern of cold tumor–related TIP signature to that of hot tumor.

In our screen, we treated MDA-MB-231 cells, a human TNBC cell line, with a library of 8199 compounds that included Food and Drug Administration (FDA)–approved drugs, epigenetic inhibitors, natural products, etc. The activity of each compound was scored according to its ability to alter the expression of our selected TIP signature genes [6 of 15 of these genes expressed in MDA-MB-231 cell lines, and the screening dataset has been described in our previous work (19)].

The highest scored compounds are ENMD-2076, JQ1, ICG001, epinephrine bitartrate, and flunixin meglumin (Fig. 4, A and B). Notably, JQ1 and ICG001 were reported to be associated with tumor immunotherapy. Multiple publications report that JQ1, a bromodomain inhibitor, could promote T cell–mediated antitumor immunity (3638). The β-catenin inhibitor ICG-001 could increase the proinflammatory signature cytokines and decreased the anti-inflammatory signature molecule interleukin-10 (IL-10) in dendritic cells and thus may change the immune phenotypes of tumor microenvironment (39).

Fig. 4 Aurora kinase inhibitors up-regulate expression of hot tumor–associated genes and down-regulate expression of cold tumor–associated genes.

(A) Potential hot tumor–boosting compounds ranked by the scores of reversing activities for the TIP signature genes. (B) The highest scored compounds are ENMD-2076, JQ1, ICG001, epinephrine bitartrate, and flunixin meglumine. (C) Analysis of chemokines mRNA expression by qPCR. Data are representative of three replicates. (D) GSEA of RNA-seq results of ENMD-2076– and TAK-901–treated MDA-MB-231 cells. NES, normalized enrichment score. FDR, false discovery rate. (E) Selected pathways and corresponding representative genes perturbed by ENMD-2076 and TAK-901 in MDA-MB-231 and 4T1 cell lines.

We picked up ENMD-2076, which scored the highest, for further study (Fig. 4, A and B). Since ENMD-2076 is an AKI, another known AKI, TAK-901, was also examined for the purpose of double checks. We found that TAK-901 also exhibited the similar pattern of gene expression as ENMD-2076 (Fig. 4C). We conducted a follow-up RNA sequencing (RNA-seq) study to examine transcriptome-wide effects of ENMD-2076 with both MDA-MB-231 cells and 4T1 cells (a murine TNBC cell line). We used gene set enrichment analysis (GSEA) and heatmap illustration to characterize the global sets of up-/down-regulated genes following ENMD-2076 treatment of each cell line (Fig. 4, D and E). As expected, ENMD-2076 could reprogram the expression of TIP signature genes (Fig. 4, D and E). Given its known ability to inhibit aurora kinases, we also found that ENMD-2076 treatment significantly inhibited AURK signaling pathway activity (Fig. 4E and fig. S2).

Furthermore, we also found that ENMD-2076 treatment inhibited the known protumor immune signal transducer and activator of transcription (STAT) pathway and promoted the activity of the antitumor immune STAT3 and T cell activation pathways as well as other immune pathways (chemokine-mediated, immune cell–associated, and inflammatory response–associated pathways; Fig. 4E and fig. S2). In light of these results, we tested another AKI, TAK-901, and RNA-seq analysis showed similar trends as with ENMD-2076 (Fig. 4, D and E). Our HTS2 screening and subsequent confirmatory experiments thus revealed that these AKIs can reprogram the expression of our selected TIP signature genes to some degree and may therefore hold potential as agents to convert immune cold tumors into immune hot tumors.

Aurora kinase A mediates the inhibitor-induced activation of TH1-type chemokine expression

To understand how AKIs regulate the expression of these TIP genes, TH1-type chemokine genes CXCL10 and CXCL11 were selected for further study for three reasons. First, our HTS2 screening data for cells treated with ENMD-2076 or quantitative polymerase chain reaction (qPCR) results of TAK-901 showed that, among our selected TIP signature genes, CXCL10 and CXCL11 were the most strongly affected by treatment (Fig. 4C). Second, our initial analysis of breast tumor dataset [METABRIC (24)] showed that high CXCL10 or CXCL11 expression is predictive for better prognosis (Fig. 3G). Last, the activation of these TH1-type chemokines has been reported to enhance the trafficking of effective T cells into tumor microenvironments and thus may enhance the IBC response in colon and ovarian cancers (6, 7). We therefore focused our attention on the direct effects of ENMD-2076 and TAK-901 on the expression of CXCL10 and CXCL11.

Confirming our results from our earlier RNA-seq analysis (Fig. 4E), we conducted additional assays, which again validated that both compounds do increase the expression of both CXCL10 and CXCL11 (Fig. 5, A and B). Moreover, these assays showed that these increases of mRNA level are dose dependent (Fig. 5, A and B). At the protein level, enzyme-linked immunosorbent assay of cell exudates showed that ENMD-2076 treatment caused a significant increase in secretion of both the CXCL10 and CXCL11 chemokines from MDA-MB-231 cells (Fig. 5C).

Fig. 5 AKIs activate TH1-type chemokine expression through inhibiting AURKA-STAT3 signaling pathway.

(A and B) Effects of ENMD-2076 (A) and TAK-901 (B) on the mRNA expression of CXCL10 and CXCL11 in breast cancer cells. (C) Protein level of CXCL10 and CXCL11 measured in MDA-MB-231 cells by enzyme-linked immunosorbent assay after ENMD-2076 treatment. DMSO, dimethyl sulfoxide. (D) The knockdown efficiency of AURKA (left) and the expression of CXCL10 and CXCL11 were quantified by qPCR (right two panels). (E) The overexpression of AURKA was confirmed by Western blot (left), and the expression of CXCL10 and CXCL11 was quantified by qPCR (right two panels). GAPDH, glyceraldehyde-3-phosphate dehydrogenase. (F) The knockdown efficiency of STAT3 (left) and the expression of CXCL10 and CXCL11 (right two panels) were quantified by qPCR. (G to I) MDA-MB-231 cells were treated with ENMD-2076 (5 μM) or interferon-γ (IFN-γ) for 24 hours, and STAT protein levels were monitored by Western blot with the indicated antibodies. (J and K) MDA-MB-231 cells were transfected with the shRNAs (J) or overexpression lentiviral vectors (K) of AURKA, and STAT3 protein levels were monitored. shNC, negative control. (l) Illustration of AKIs enhancing CXCL10/CXCL11 expression through inhibition of STAT3 phosphorylation. *P < 0.05,**P < 0.01, and ***P < 0.001 (Mann-Whitney U test). Error bars depict SEM.

Given that both ENMD-2076 and TAK-901 are known AKIs, we next sought to determine whether aurora kinase enzymes have a direct functional role in the observed inhibitor-induced TH1-type chemokine activation. We therefore used a lentivirus-based method to introduce a short hairpin RNA (shRNA) that successfully knocked down the expression of aurora kinase A (AURKA) and found that AURKA knockdown significantly increased the expression of CXCL10 and CXCL11 (Fig. 5D). We also overexpressed AURKA in MDA-MB-231 cells, which caused a significant decrease in the expression of CXCL10 and CXCL11 (Fig. 5E).

In addition, investigation of genomic data from patients with breast cancer showed that ARUKA was frequently genetically amplified in most of patients (fig. S3A; the patient-ranked list in this figure is the same with Fig. 3C), which conversely exhibited relatively low expression and frequent deletion of CXCL10 and CXCL11 (Fig. 3C). Moreover, risk analysis indicated strong AURKA expression was common in the patients with breast cancer with high-risk scores and short OS times (fig. S3B; the upper and middle panels in this figure are the same with Fig. 3H). Kaplan-Meier analysis further indicated that high AURKA expression in tumors was significantly correlated with poor prognosis in patients with breast cancer (fig. S3C). Moreover, we observed that the expression levels of CXCL10 and CXCL11 were inversely correlated with AURKA in patients with breast cancer (fig. S3D), especially in patients with TNBC (fig. S3, E and F).

Together, our experimental results demonstrate that AURKA does function in ENMD-2076– and TAK-901–induced activation of the TH1-type chemokines CXCL10 and CXCL11, and our clinical data analysis establishes a link between increased AURKA levels, decreased CXCL10 and CXCL11 levels, and worse clinical outcomes in patients with breast cancer.

STAT3 phosphorylation mediates the AURKA inhibition–induced activation of TH1-type chemokine expression

To further determine the downstream targets of AURKA in the regulation of TH1-type chemokine expression, we examined these significantly altered signaling pathways in the cells treated by these two inhibitors. Our RNA-seq analysis of breast cancer cells treated with ENMD-2076 and TAK-901 showed that these compounds significantly induce the expression of multiple genes in STAT3 signaling pathway (Fig. 4E), so we investigated whether STAT3 may also function in the inhibitor-induced activation of CXCL10 and CXCL11. A GSEA of our RNA-seq data conformed significant enrichment in the inhibitor-treated samples for known STAT3 target genes (fig. S4A) (40). Moreover, a clinical outcome analysis indicated that strong STAT3 expression in breast tumors is correlated with poor prognosis (fig. S4B).

To experimentally determine whether the inhibitor-induced activation of CXCL10 and CXCL11 expression is mediated by STAT3, we used two shRNAs to knock down STAT3 expression in MDA-MB-231 cells. Knockdown of STAT3 resulted in significantly increased CXCL10 and CXCL11 expression (Fig. 5F). We also conducted assays, which showed that ENMD-2076 treatment clearly inhibited the phosphorylation of STAT3, which is known to be required for STAT3 functional activation (Fig. 5G) (41, 42). These effects are apparently STAT3 specific, as ENMD-2076 treatment showed almost no effect on the expression levels or phosphorylation status of other STAT family members such as STAT1, STAT5, or STAT6 (Fig. 5, H and I).

To assess whether the function of STAT3 is regulated by AURKA, the STAT3 phosphorylation is examined in TNBC cells. We used two shRNAs to knock down AURKA expression in MDA-MB-231 cells, which clearly leads to the decreased phosphorylation of STAT3 (Fig. 5J). In addition, the AURKA protein is overexpressed by lentivirus and the STAT3 phosphorylation is significantly increased (Fig. 5K). These analyses and experiments together establish that the alteration of the phosphorylation status of STAT3 is regulated by AUKRA and involved in the aurora inhibitor–induced activation of TH1-type chemokines in TNBC cells (Fig. 5L).

Aurora inhibitors elevate TH1-type chemokine levels and promote CD8+ T cells infiltration in vivo

Since the CD8+ T cell infiltration into the tumor microenvironment is a key biomarker to define cold and hot tumors, we checked whether inhibitor-induced expression of CXCL10 and CXCL11 in breast cancer cells could enhance the infiltration of CD8+ T cells. We found that ENMD-2076 treatment of MDA-MB-231 cells significantly increased the migration of in vitro activated CD8+ T lymphocytes into the growth medium (Fig. 6, A and B); a similar increase was also observed following the positive control interferon-γ (IFN-γ) treatment (Fig. 6, A and B). We further conducted an experiment that showed, consistent with the known role of CXCL10 and CXCL11 in CD8+ T cell recruitment, that this aurora inhibitor–induced migration of T cells was blocked when the cells were pretreated with an antibody known to disrupt the ability of the CXCR3 receptor to sense the CXCL10 and CXCL11 chemokines (Fig. 6, A and B). Thus, these results demonstrated that the AKIs enhance the infiltration of CD8+ T cells by activating the expression of the TH1-type chemokines CXCL10 and CXCL11 in TNBC cells (Fig. 6C).

Fig. 6 AKIs improve the infiltration of CD8+ T cell and the efficacy of anti–PD-1 immunotherapy in TNBC in vivo.

(A) Representative flow cytometric analysis of CD8+ T cells migration in vitro. (B) The percentage of CD8+ T cell migration described in (A) from independent sample (Mann-Whitney U test). (C) Illustration of AKIs enhancing CD8+ T cell chemotaxis. (D) Illustration of animal models. 4T1 cells were injected into the mammary fat pads of mice. Animals were administrated when the volumes of tumors were about 100 mm3. (E) Representative images of Immunohistochemistry (IHC) staining showing the distribution of CD8+ T cells. Scale bar, 50 μm. P, E, and T represent PD-1–, ENMD-2076–, and TAK-901–treated groups, respectively; P+E and P+T represent particular combination therapies, respectively. (F) Representative IHC analysis of CD8+ T cells in independent samples (Mann-Whitney U test). (G) The expression level of Cxcl10 and Cxcl11 in 4T1 tumors was analyzed by qPCR; each square represents data from an independent mouse. (H) The growth of 4T1 tumors was measured by tumor volume; volume (mm3) = [width2 (mm2) × length (mm)]/2. (I) Tumor sizes at day 25 (sample-paired Student’s t test). *P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001. Error bars depict SEM.

To evaluate the effect of ENMD-2076 and TAK-901 treatment on the infiltration of CD8+ T cells and cancer immunotherapy against TNBC in vivo, we injected murine 4T1 cells into mammary pads to generate a syngeneic murine model. Once the size of these TNBC tumors reached about 100 mm3, the mice were randomly divided into six groups for different treatments: anti–PD-1 antibody alone, ENMD-2076 alone, TAK-901 alone, combination anti–PD-1 + ENMD-2076, combination anti–PD-1 + TAK-901, and an immunoglobulin G (IgG) control group (Fig. 6D).

We used immunohistochemistry staining of tumor sections to examine the distribution and trafficking of CD8+ T cells. As expected, AKIs caused increases in the infiltration of CD8+ T cells to tumor sites in contrast to control (Fig. 6, E and F). Whereas the anti–PD-1 treatment alone almost had no effects on CD8+ T infiltration, the combination therapies, anti–PD-1 + ENMD-2076 and anti–PD-1 + TAK-901, both significantly enhanced the trafficking of CD8+ T cells into the tumor microenvironment (Fig. 6, E and F). Moreover, flow cytometry analysis of fresh tumors confirmed that the proportions of intratumoral CD8+ T cells increased significantly in response to the combination of anti–PD-1 and AKI treatments, and the single compound treatment also had the increasing trend (fig. S5A). Consistently, the in vivo results also showed that ENMD-2076 treatment alone or in combination with anti–PD-1 treatment increased Cxcl10 and Cxcl11 expression in TNBC tumors (Fig. 6G).

As IFN-γ plays an important role during antitumor immune response, a previous study also proved that CD8+ T cell–secreted IFN-γ and its associated dendritic cell–derived IL-12 sensitized tumors to anti–PD-1 therapy (43), we next examined the levels of IFN-γ in tumor-infiltrating CD8+ T cells. Flow cytometry analysis of tumors revealed a significantly increased production of IFN-γ in the combination therapy group relative to anti–PD-1 alone, while the treatment of TAK-901 alone also showed an increasing trend (fig. S5B).

In addition to examine CD8+ T cells in tumors, we also checked the effect of these AKIs on immune cells in mouse spleens. By means of flow cytometry analysis, we observed higher percentages of CD8+ T cells in spleens from mice under the treatment of ENMD-2076 alone or combination of anti–PD-1 + ENMD-2076 compared with treatment of IgG alone (fig. S6A), while the other immune cell populations, including CD4+ T cells, macrophages, and myeloid-derived suppressor cells, exhibited no significant differences (fig. S6, B to D). Collectively, these results indicated that the AKIs increase the expression of CXCL10 and CXCL11, enhance the infiltration of effective CD8+ T cells into the tumor microenvironment, and partially reprogram the cold immune microenvironment of TNBC to hot, while showing no inhibition effect on immune cell populations in spleen.

Aurora inhibitors enhance the efficacy of ICB therapy on TNBC in vivo

To evaluate the effect of ENMD-2076 and TAK-901 treatment on the ICB response in TNBC, tumor size of these six groups of mice was measured every 5 days throughout the course of a 25-day drug treatment (Fig. 6D). As would be expected as a cold tumor, and consistent with previously reported basic and clinical observations of TNBC tumors, we found that treatment with anti–PD-1 alone had almost no effect on TNBC tumor growth in our model (Fig. 6, H and I).

Since AKIs are reported as cell mitotic inhibitors and thus inhibit cell proliferation (44), we found that treatment with ENMD-2076 or TAK-901 alone did inhibit tumor growth as expected (Fig. 6, H and I). However, the combination treatment of anti–PD-1 + ENMD-2076 exerted a significantly more pronounced tumor growth inhibition effect than ENMD-2076 alone or anti–PD-1 alone. Consistently, the combination treatment of anti–PD-1 + TAK-901 showed a similar trend in inhibiting TNBC tumor growth than anti–PD-1 alone or TAK-901 alone (Fig. 6, H and I). Thus, our results clearly suggested that the combination of these AKIs significantly enhances the efficacy of anti–PD-1 treatment on TNBC tumor growth in vivo.

DISCUSSION

Here, we provide a new approach to effectively identify enhancing agents for immunotherapy, which combines a TIP gene signature and HTS2. This screen and subsequent experiments identified that treatment with two AKI drugs, ENMD-2076 and TAK-901, up-regulates the expression of a number of antitumor chemokines in both human and murine TNBC cells, promote the infiltration of effective T cells into the tumor microenvironment, and thus significantly sensitized TNBC cells to anti–PD-1–based immune checkpoint therapy in preclinical models. Furthermore, we show that the TH1-type chemokines CXCL10 and CXCL11 are strongly induced by ENMD-2076 and TAK-901 and provide both genetic and biochemical evidence that both AURKA and STAT3 phosphorylation are involved in this inhibitor-induced activation of TH1-type chemokine expression (Fig. 7).

Fig. 7 Schematic illustration of the novel approach for the discovery of combination immunotherapy agents and the mechanism of AKIs enhancing immunotherapy efficacy.

It is quite challenging for the potential compounds, which could be used in the combination immunotherapy, by traditional high-throughput drug screening. There have not been any high-throughput screens reported to date in this area. This is due, in part, to the lack of appropriate in vitro cell models: consider that such a screen would, rather than simply examining postcompound-treatment cell viability, need to monitor the complex process of immune cell infiltration. As a strategy to overcome these challenges, we used a drug screening method named HTS2 that uses gene expression signatures rather than cellular phenotypes as readouts. To facilitate HTS2 screening, we also successfully identified a set of TIP genes, which could be used to distinguish hot with cold tumors, as readouts. We believe that our method of combining TIP gene signature and HTS2 might be a promising high-throughput approach for the discovery of immunotherapy-enhancing agents and has great potential to promote cancer immunotherapy.

ENMD-2076 is an orally bioavailable small-molecule inhibitor of AURKA as well as other kinases such as vascular endothelial growth factor receptors and fibroblast growth factor receptors. This drug is active against preclinical TNBC models, but a phase 2 clinical trial reported that, when used alone, it elicits only partial therapeutic responses and limited clinical benefit in patients with TNBC (45). In addition, this compound has been tested as a single agent in other phase 2 clinical trials against ovarian clear cell carcinoma (46) but likely did not meet the preset bar for efficacy. TAK-901 is a multitargeted AKI. In preclinical models, TAK-901 exhibited potent activity against multiple human solid tumor types; TAK-901 has been tested in two phase 1 clinical trials, and we did not find any record of ongoing or planned phase 2 clinical trials for this drug. Viewed collectively, it thus looks like the therapeutic efficacy of these drugs as single agents is not so promising. Notably, our study illustrates an alternative clinical practice opportunity for both ENMD-2076 and TAK-901: as the small-molecule component alongside checkpoint blockade antibodies in combination immunotherapies.

Our investigations of the molecular signaling pathways associated with ENMD-2076 and TAK-901 treatment revealed that the activation of TH1-type chemokine genes CXCL10 and CXCL11 involves AURKA and STAT3. These discoveries are consistent with previous studies, which reported that AURKA inhibition can significantly decrease STAT3 phosphorylation and inhibit its function (47) and which reported that activated (phosphorylated) STAT3 down-regulates the expression of TH1-type chemokines (48). It also bears mention that AURKA is a frequently amplified and overexpressed gene in multiple cancers, including TNBC, and STAT3 is known to promote proinflammatory signaling (49, 50). STAT3 has been shown to inhibit nuclear factor κB binding to TH1-type genes in growing tumors, including in tumor-infiltrating immune cells (51). In addition to aurora kinase, our results indicate that STAT3 might be another promising target for the combinational immunotherapy and need to be studied further. Future studies should also focus on elucidating the detailed mechanisms that connect AURKA, STAT3 phosphorylation, and the expression of TH1-type chemokines. A clear understanding of these signaling interactions should enable the design of new therapeutic strategies that exploit the highly impactful role of immune cell infiltration on the efficacy of immunotherapy.

Our study provides a new approach by combining TIP gene signature and HTS2 to effectively identify enhancing agents for immunotherapy. Applying this new approach, we identify and demonstrate that the AKIs ENMD-2076 and TAK-901 can enhance the efficacy of checkpoint blockade immunotherapy in preclinical TNBC models and points the way forward for a new combination immunotherapy strategy that might improve the response rates of the unfortunate patients suffering from the frequently immunotherapy-resistant TNBC, a highly aggressive disease with extremely poor prognosis and few current treatment options.

MATERIALS AND METHODS

Study design

We identified the top 15 TIP signature genes using IE analysis and used our HTS2 approach (18, 19) to screen 8199 molecules by monitoring their effects on the expression of a defined set of genes. RNA-seq and GSEA analysis were used for investigation of candidate compounds. All the in vitro and in vivo studies were performed using the MDA-MB-231 cell line, the 4T1 cell line, or Balb/c mice to detect candidate compounds that promote effective T cell infiltration into the tumor microenvironment and that improve anti–PD-1 efficacy in controlling tumor growth.

TIP signature

To obtain a set of genes that could represent the tumor immune phenotype, we screened all immune and immune therapy–related genes from 30,000 literatures using IE, a text-mining approach (22). In total, 15 genes were ranked at the top of the immune and immune therapy–related gene list and were subsequently classified into either the hot tumor– or cold tumor–related category on the basis of the classification from the published review (4). Information of patients with melanoma (n = 266) was collected from this paper (23) that has already classified a tumor into either a hot or cold tumor. According to the review (4), cold tumor–related genes (such as CXCL1) should be lowly expressed in hot tumor; then, if they were lowly expressed in the test melanoma cohorts (23), a +1 score was assigned to them, vice versa, a −1 score was assigned. Similarly, hot tumor–related genes (such as CXCL10) should be highly expressed in hot tumors; if they were highly expressed in the test melanoma cohorts (23), a +1 score was assigned to them, vice versa, a −1 score was assigned. Individual patients that own these 15 genes could obtain the sum of these scores. Hence, patients with scores above 0 were predicted as hot tumor–related patients and vice versa. Next, the predicted hot or cold tumors were pairwise compared to the actual hot or cold tumors (23); the prediction accuracy and ROC (receiver operating characteristic) curve could be obtained according to the comparison results. The R package of “pROC” was used to generate an ROC curve and an AUC (area under curve) value. To check the expression level of each gene in each patient, the unsupervised clustering was only performed on genes (rows) rather than samples (columns). The correlation analysis for the infiltration level of CD8+ T cells and the expression level of six TIP signature genes were performed with tumor immune estimation resource (TIMER) (15).

Clinical data processing

The breast cancer datasets were collected from the “cBioPortal for Cancer Genomics” website [www.cbioportal.org (24)]. Pearson correlation analysis between different genes was performed on the basis of the expression level of these genes in 1980 patients with breast cancer. The Kaplan–Meier (KM) plots for genes including CXCL10, CXCL11, CCL5, CXCR3, CXCL1, CXCL2, AURKA, and STAT3 were plotted on the website of Kaplan-Meier Plotter (http://kmplot.com/analysis/). Analysis of the mutation of individual genes (CXCL9, CXCL10, CXCL11, CXCR3, CXCR4, CCL5, CD8A, CD8B, CD4, CD3E, CD274, PDCD1, CXCL1, CXCL2, and AURKA) and the gene signature were performed on the “cBioPortal for Cancer Genomics” website on the basis of the mutation information of patients with breast cancer [METABRIC (24)].

Performance evaluation of TIP and other published biomarkers on predicting ICB response

To evaluate the performance of TIP-associated genes on predicting ICB response, RNA-seq datasets from four cohorts were obtained from the original paper (2629). For each dataset, we standardized the transcriptome data across patients by quantile normalization and further normalized the expression values of each gene by subtracting the average among all samples (z score). Next, the z score was used as the input of downstream visualization: Bar plot was used to represent the averaged z score of hot tumor–related genes, and heatmap was used to visualize the clustered results of each patient in these four datasets with this z score by an R package of GENE-E. We also compare the response prediction of TIP with other biomarkers published in the literature by following the methods provided by the updated version of TIDE (16, 30) over seven datasets (27, 28, 3135).

Clinical outcome analysis

To assess the risk score of each patient, the breast cancer dataset [METABRIC (24)] was used. First, the value of coefficient for each gene was computed on the basis of the expression levels and the living months as well as the living status (dead or not) of each gene in 1980 patients with “coxph test” from the “survival” R package. Second, for individual patient, a risk score was obtained by sum of the value of coefficient*expression level of each gene.

To explore the clinical relevance of the expression of TIP signature genes with patients with cancer over 28 TCGA cancer types, the multivariable Cox proportional hazard model was performed with the “survival” module of TIMER website (15). Briefly, the Cox’s regression model was defined with the formula of Surv(CancerType) ~ variables (only hot tumor–related genes were included here). The model was fitted with the function of coxph() from R package “survival,” and the coefficient “coef” was used as a regression coefficient.

Drug screening and enrich score calculation

The drug-screening data used in this project were generated by Shao et al. (19). The method of calculating the activity score to reverse cold tumor into hot tumor for each compound was similar with the method described in the original connectivity map literature (17) and in our previous work (19).

Six of 15 TIP genes (those expressed in MDA-MB-231 cell lines and with available screening probes) including CXCL10, CXCL11, CXCL1, CXCL2, CCL20, and CXCR4 were combined to serve as the “query” gene set, the HTS2 expression profiles of all screened compounds served as the “target” gene sets, and Kolmogorov-Smirnov test was used to compare the similarity of each paired query and target to assign a final similarity score for each compound.

In detail, the R script were downloaded from the GitHub of the original paper (https://rdrr.io/bioc/gCMAP/src/R/connectivity_score-methods.R) (17) and modified with R. All screened drugs were then ranked according to their ability to alter the expression of those six genes. For each drug, compute an enrichment score [(ksi up) and (ksi down)] for the normalized counts of probes representing the up-regulated (CXCL1, CXCL2, and CCL20, need to be down-regulated after drug treatment) or down-regulated (CXCL10, CXCL11, and CXCR4, need to be up-regulated after drug treatment) genes in the signature. Set n as the total number of genes that included for screening (1224) and t (6, here) as the number of genes in the signature. Then, a vector (V) of the position (1 to n) of each gene in the ordered list of all genes was constructed, followed by computing two valuesa=maxj=1t[jtV(j)n]b=maxj=1t[V(j)n(j1)t]and if a > b, set ksi = a, or if b > a, set ksi = −b. The enrich score was set to zero when (ksi up) and (ksi down) have the same algebraic direction. Otherwise, set si = (ksi up) – (ksidown), set p = max(si) and q = min(si) across all drugs. The enrich score Si is defined as si/p where si > 0, or –(si/q) where si < 0. All drugs were then ranked in descending order of S and ksup (17).

RNA sequencing

A total of 2 × 105 MDA-MB-231 cells were plated in a six-well plate for 24 hours and then treated with ENMD-2076 at 2.5 or 5 μM (or TAK-901 at 1 or 2 μM) for 24 hours. A total of 1 × 105 4T1 cells were plated in a six-well plate for 24 hours and then treated with ENMD-2076 at 1 or 2 μM (or TAK-901 at 1 or 2 μM) for 24 hours. Cells were harvested, and total RNA was isolated using TRIzol (Invitrogen). The libraries were constructed by using the NEBNext Ultra RNA Library Prep Kit for Illumina (New England Biolabs) and sequenced on a HiSeq 2500 sequencing system (Illumina).

RNA-seq FASTQ data were mapped to the reference genome (hg38 and mm10) using Bowtie/TopHat. The reads were counted with HTSeq (52), and the differential expression between experimental groups was quantified using edgeR (53) (fold change > 2 and P < 0.05) and DESeq2 (54).

Gene Ontology and GSEA

Gene Ontology (GO) analysis was performed in DAVID v6.7 (https://david.ncifcrf.gov/) and Gorilla (http://cbl-gorilla.cs.technion.ac.il) with differentially expressed genes (DEGs) (fold change > 2 and P <0.05). GSEA was performed using the GSEA v3.0 software. First, the complete tables of genes ranked according to their fold change [drug versus dimethyl sulfoxide (DMSO)] were served as the input datasets. Second, four different gene sets that were obtained from the studies of (4) and (55) served as input gene sets, the four gene sets containing (i) the up-regulated (CXCL10, CXCL11, CXCL13, CXCL14, CXCL16, CXCL17, CXCL6, CXCL9, CXCR3, CXCR5, CCL20, CCL24, TNF, and CCR6) or down-regulated genes (CXCL1, CXCL12, CXCL2, CXCL3, CXCL5, CXCL8, CXCR1, CXCR2, CXCR4, CXCR6, CCL1, CCL11, CCL13, CCL17, CCL2, CCL25, CCL26, CCL27, CCL28, and CCL7) in hot tumor (4), (ii) up-regulated genes with STAT3 knocked down, (iii) STAT3 target genes detected by chromatin immunoprecipitation sequencing, and (iv) overlapped genes of (i) and (iii) (55). Third, the enrichment score and false discovery rate value were applied to sort each dataset after gene set permutations were performed 1000 times for the analysis.

Plasmids and antibodies

Lentiviral shRNAs were provided by the Vector Core at Tsinghua University. The scrambled shRNA used was GGTGTATGGGCTACTATAGAA. To construct the human AURKA ectopic expression plasmid, plvx-AURKA (a full-length human AURKA) was cloned from MDA-MB-231 complementary DNA (cDNA). Antibodies used included monoclonal anti-human AURKA (1:2000, Cell Signaling Technology), anti-human STAT1 (1:2000, Cell Signaling Technology), anti-human p-STAT1-Tyr701 (1:2000, Cell Signaling Technology), anti-human p-STAT1-Ser727 (1:2000, Cell Signaling Technology), anti-human STAT5 (1:2000, Cell Signaling Technology), anti-human STAT6 (1:2000, Cell Signaling Technology), anti-human STAT3 antibody (1:2000, Cell Signaling Technology), anti-human p-STAT3 (1:2000, Cell Signaling Technology), anti-human glyceraldehyde-3-phosphate dehydrogenase (GAPDH) antibody (1:2000, Cell Signaling Technology), anti-mouse CD8 (1:2000, Abcam), anti-human CXCR3 blocking antibody (BioLegend), anti-mouse CD4 Super Bright 600 (eBioscience), anti-mouse CD8a PE-eFluor 610 (eBioscience), anti-mouse CD3e PerCP-Cyanine 5.5 (eBioscience), anti-mouse IFN-γ allophycocyanin (APC) (eBioscience), anti-mouse CD45 eFluor 450 (eBioscience), anti-mouse Ly-6G/Ly-6C Super Bright 780 (eBioscience), anti-mouse F4/80 PE-CY7 (eBioscience), anti-mouse CD11b fluorescein isothiocyanate (FITC) (eBioscience), anti-human CD4 PE (eBioscience), anti-human CD8 FITC (eBioscience), and anti-human CD3e APC (eBioscience).

Cell culture

The 4T1 and MDA-MB-231 cell lines were obtained from the Institute Cancer Cell Line Encyclopedia of China. 4T1 cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM) (Gibco) supplemented with 10% fetal bovine serum (Gemini) and 1% penicillin/streptomycin (Gibco). MDA-MB-231 cells were cultured in R1640 (Gibco) supplemented with 10% fetal bovine serum (Gemini) and 1% penicillin/streptomycin (Gibco). Cells were cultured at 37°C with a 5% CO2 atmosphere. The cell lines were tested for mycoplasma contamination.

Chemicals

The screening library (1154 FDA-approved drugs) as well as ENMD-2076 (S1181) and TAK-901 (S2718) were purchased from Selleck. Other chemicals were supplied by Shanghai Institute of Materia Medica. All chemicals were dissolved in DMSO for the in vitro studies.

Mice and reagents

The laboratory animal facility has been accredited by AAALAC (Association for Assessment and Accreditation of Laboratory Animal Care International) and the IACUC (Institutional Animal Care and Use Committee) of Tsinghua University approved all animal protocols used in this study.

BALB/c female mice (5 to 6 weeks old) were housed in the animal facility of the Laboratory Animal Research Center of Tsinghua University under specific pathogen–free conditions. All the animal experiments were approved by Tsinghua University’s Animal Care and Use Committee. Briefly, subconfluent 4T1 cells (5 × 105) were suspended in 0.04-ml DMEM and injected into the mammary fat pads of mice under anesthesia with avertin. Animals were administrated when the volumes of tumors were about 100 mm3. Mice were treated with ENMD-2076 at 100 mg/kg orally twice weekly for 3 weeks, with IgG, or with anti-mouse PD-1 at 10 mg/kg intraperitoneally (ip) twice weekly for 3 weeks. Mice were treated with TAK-901 treatment at 30 mg/kg ip twice weekly for 3 weeks. Tumor volumes were measured with digital calipers and were calculated with the formula: volume (mm3) = [width2 (mm2) × length (mm)]/2. Tumor growth was detected by bioluminescence imaging with an IVIS Lumina II instrument (PerkinElmer). After 25 days from the beginning of treatment, mice were euthanized. Tumors were then isolated for fluorescence-activated cell sorting (FACS). Part of each tumor was fixed in 4% paraformaldehyde for immunohistochemistry.

Quantitative polymerase chain reaction

RNA that was exacted with TRIzol (Invitrogen) was used for cDNA synthesis with the Thermo Scientific RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific). qPCR was performed using a StepOne Plus real-time PCR system (KAPA). The specific primers used for qPCR are listed in table S1.

T cell migration assays

In vitro migration assays were performed in a Transwell system (Corning) with a polycarbonate membrane of 6.5-mm diameter with a 3-mm pore size as described previously (56). MDA-MB-231 cells were seeded in the lower chamber of the Transwell insert with ENMD-2076 treatment (5 μM) for 48 hours. Then, activated T cells pretreated with antihuman CXCR3, or isotype were added to the upper chamber. The phenotype and number of T cells in the lower chambers were measured with FACS (LSR II, BD Biosciences) after incubated at 37°C for 6 hours.

Immunoblotting

Cells were lysed with radioimmunoprecipitation assay lysis buffer containing phenylmethylsulfonyl fluoride (1 mM) and 1× protease inhibitor cocktail (Biotool, Beijing). After bicinchoninic acid (BCA) protein quantification (57), protein samples were separated by SDS–polyacrylamide gel electrophoresis and transferred to polyvinylidene membranes (Bio-Rad) and probed with antibodies indicated in the legends of each relevant figure. The bands were visualized by using SuperSignal West Pico reagent (Thermo Fisher Scientific). GAPDH was used as an internal control.

Statistical analysis

Two-tailed unpaired Student’s t tests were used to compare the differences of two groups; log-rank tests were used for Kaplan-Meier survival analysis; Pearson’s coefficient correlation tests were used for gene correlation analysis; chi-square test was used for association analysis between patients with inversely correlated expression of AURKA and CXCL10/ CXCL11 and patients with TNBC. The level of significance was set at *P < 0.05, **P < 0.01, and ***P < 0.001.

Code availability

All codes used in this study followed the manuals of each R package, which are fully open accessed.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/7/4/eabd7851/DC1

https://creativecommons.org/licenses/by-nc/4.0/

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

References and Notes

Acknowledgments: We thank the staff members of the Animal Facility of the Laboratory Animal Research Center of Tsinghua University. Funding: This work was financially supported by National Natural Science Foundation of China (grant numbers 31660254, 81673460, and 31371314). Author contributions: D.W., Haiyan Wang, and S.L. jointly conceived the idea of the experiment and conception. Haiyan Wang, Q.W., and Z.J. prepared samples and conducted measurements. S.L. analyzed all datasets and processed sequencing data. All authors coordinated the work. D.W., S.L., Haiyan Wang, and Q.W. wrote the manuscript with assistance from all the other authors. Competing interests: Haiyan Wang, S.L., and D.W. are inventors on a patent application related this work filed by Tsinghua University (filed: 31 December 2018; no. 201811651274.1). The authors declare that they have no other competing interests. Data and materials availability: All raw and processed RNA-seq data for ENMD-2076– and TAK-901–treated cell lines were deposited at the GEO database under accession number GSE160071. All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article