Single-cell RNA-sequencing on mammary epithelial cells across four developmental stages (nulliparous, mid-gestation, lactation, and post involution) were described by the original report (Bach et al., 2017). We accessed the processed data of this scRNA-seq dataset from the Gene Expression Omnibus under accession number GSE106273 and raw data of the scRNA-seq from Sequence Read Archive (SRA) under accession number PRJNA416110 for our analysis. scRNA-seq data on different mouse models of breast cancer were retrieved from GEO under accession number GSE123366 and SRA under accession number PRJNA508501 (Yeo et al., 2020). scRNA-seq data on primary TNBC were retrieved from GEO under accession number GSE118389 and SRA under accession number PRJNA485429 (Karaayvaz et al., 2018), and scRNA-seq data on brain metastatic TNBC were retrieved from GEO under accession number GSE143423 and SRA under accession number PRJNA600483.

RStudio software package Seurat was used for the analysis of matrix data of scRNA-seq. Sample objects were integrated or merged by either sample individuals or their disease status. Anchors were generated, listing all the Seurat objects as input, and integrated/merged data were created for integration. In the integrated/merged Seurat object, “nFeature_RNA,” “nCount_RNA,” and “percent.MT” were defined to represent number of genes, number of transcripts, and percentage of mitochondrial initiated genes in each cell of each sample, respectively. Quality control (QC) was performed by removing the low-quality events and doublets by excluding the cells with extreme nCount_RNA or nFeature_RNA and high percent.MT.

Canonical cell-type markers were used to identify and define the major cell types. For the epithelial/cancer cell lineage, we used common epithelial cell marker genes [CDH1 (Cdh1), KRT18 (Krt18), and KRT5 (Krt5)] and endothelial cell marker genes [CDH5 (Cdh5), PECAM1 (Pecam1), and EMCN (Emcn)]. Total immune cells were identified by PTPRC (Ptprc), CD68 (Cd68), and CSF1R (Csf1r). For the mesenchymal cell component, COL1A1 (Col1a1), COL1A2 (Col1a2), COL3A1 (Col3a1), and FAP were used. The cancer/epithelial cell clusters were extracted from other cell types and re-clustered for further analysis.

The scRNA-seq data of epithelial/cancer cells were loaded into Slingshot packages, and major clusters were identified for pseudotime analysis. Analysis on inferred trajectory and principal curves demonstrated smoothed representations of each subcluster. Output principal curves were smoothed representations of each disease status, and pseudotime values were computed by projecting the cells onto the principal curves. For the Ingenuity Pathway Analysis (IPA), differentially expressed genes of each disease status were identified and were loaded into the IPA client (QIAGEN) to demonstrate the key ingenuity canonical pathways and upstream regulators. The ingenuity canonical pathways, upstream regulators, and graphical summary of the regulatory networks were exported from the IPA.

MDA-MB-231 cell (RRID: CVCL_0062) was purchased from American Type Culture Collection (ATCC) and was cultured in L-15 medium supplemented with 10% FBS, penicillin (100 U/ml), and streptomycin (100 μg/ml). The cells have been authenticated using STR profiling within the last 3 years. All experiments were performed with mycoplasma-free cells. MDA-MB-231 cells were cultured in six-well plates to 30–60% confluency and shRPL27A (Santa Cruz Biotechnology), or scramble lentiviral particles were mixed in complete medium together with polybrene (10 μg/ml). Cells were cultured with lentiviral particles for 4 h, and then the medium was changed. Cells were cultured until confluency and were passaged for further use. For the EIF2 inhibitor assay, 20 and 50 μm Sal003 were used, and the same amount of DMSO was used as controls.

The RNeasy Mini kit (QIAGEN) was used to extract total RNA from samples following transfection according to the manufacturer’s protocols. The expression of human RPL27A was determined using the SYBR Green Master Mix kit (Roche). Human GAPDH was used as an internal control gene, and relative expression levels were calculated by using the 2^–ΔΔCt method. The sequences of specific primer pairs are described below: RPL27A-F 5′-ACGGGTGAATGCTGCTAAAA-3′ and RPL27A-R 5′-GAAGAATTTGGCCTTCACGA-3′; RPL37A-F 5′-ACATGGCCAAACGTACCA-3′ and RPL37A-R 5′-TGCTG GCTGATTTCAATTTT-3′; RPL30-F 5′-ATGGCCAAACGTA CCAAG-3′ and RPL30-R 5′-AAGTGTACTTGGCGTGCTG-3′; RPL39 primer (Bio-Rad, qHsaCED0038310); RPL36 primer (Bio-Rad, qHsaCED0038433); RPL36A primer (Bio-Rad, qHsaCED0044302); RPL15 primer (Bio-Rad, qHsaCED0038184); RPL21 primer (Bio-Rad, qHsaCED0038684); RPL38 primer (Bio-Rad, qHsaCID0038174); RPS7 primer (Bio-Rad, qHsaCID0038168); RPS15A primer (Bio-Rad, qHsaCED0038046); RPS17 primer (Bio-Rad, qHsaCED0041771); RPS21 primer (Bio-Rad, qHsaCID0023899); RPS27A primer (Bio-Rad, qHsaCED0038707); RPS29 primer (Bio-Rad, qHsaCED0038808); and GAPDH-F 5′-AGCCACATCGCTCAGACAC-3′ and GAPDH-R 5′-GCCCA ATACGACCAAATCC-3′.

RIPA buffer (Thermo Fisher Scientific) and Protease and Phosphatase Inhibitor (Thermo Fisher Scientific) were used to extract proteins. The primary antibodies against antibody for RPL27A (NBP2-38025, Novus Biologicals) and GAPDH (NB300-221, Novus Biologicals) and secondary antibodies were anti-mouse IgG-HRP (7076S, Cell Signaling Technology) and anti-rabbit IgG-HRP (7074S, Cell Signaling Technology). The protein blots were visualized using an ECL Western Blotting Substrate (Thermo Fisher Scientific).

Twenty-four-well chambers (Corning) were used or migration assay, and 24-well Matrigel invasion chambers (BD Biosciences) were used for invasion assay. MDA-MB-231 cells (5 × 10⁴) were placed in 100-μl serum-free medium in the migration/migration chamber, and 750 μl complete medium was placed into the lower wells. After culturing for 24 h, cells were fixed with 100% methanol for 20 min and stained with Trypan blue for 30 min. Non-migrating cells on the upper side of the filter were removed with cotton swabs. Migration and invasion were quantified by counting the number of cells on the lower surface of the filter.

The cell number of each cluster or each disease status was determined in Seurat. Student’s two-tailed t-test was used for comparing differences between two groups, and p < 0.05 was considered statistically significant: ^∗p < 0.05; ^∗∗p < 0.01; ^∗∗∗p < 0.001, ^****p < 0.0001. The cell numbers and cell percentages were calculated and were plotted using GraphPad Prism 8.

All the scRNA-seq data are available in GEO (accession numbers are above), and the codes regarding to the data analysis are accessible upon request.

To gain insights into the heterogeneity within mammary tumors, a recent study performed single-cell transcriptional profiling of breast tumor cells from four mouse models (Yeo et al., 2020). The MMTV-Neu (NEU) tumor was ER^–/PR^–/HER2⁺, while the MMTV-PyMT (FF99WT) tumor was ER^–/PR^–/HER2^low; the BRCA1-null (BRCA1) tumor which mimics basal-like breast cancers was ER^–/PR^–/HER^–, and the 4T1 transplant tumor was a known ER^–/PR^–/HER^– mouse model. To study the gene profiles of the tumor cells in these four models, we reaccessed the scRNA-seq data from GEO. We first integrated all the scRNA-seq data of different samples and models and corrected the batch effects. Cell qualities in each sample were visualized by violin plots. Cells that met the requirement were retained, and cells with extreme nFeature_RNA and nCount_RNA or high percent.MT were removed (Supplementary Figures 1A,B). After clustering, the distribution of each sample (Supplementary Figure 1C) and each model (Supplementary Figure 1D) was visualized by Uniform Manifold Approximation and Projection (UMAP). In total, 14 major cell clusters were identified, and the clusters of each model were also visualized (Supplementary Figure 1E).

We then checked the expression of the canonical cell-type markers: cancer/epithelial cell marker genes (Cdh1, Krt18, and Krt5); mesenchymal cell marker genes (Col1a1, Col1a2, and Col3a1); immune cell marker genes (Ptprc, Cd68, and Csf1r); and endothelial cell marker genes (Pecam1, Emcn, and Cdh5) (Supplementary Figure 1F). Based on the transcripts of these genes, we accordingly identified the cancer cells, mesenchymal/endothelial cells in the integrated data, and each tumor model (Supplementary Figures 1G,H). Then the cancer cells were subset from the integrated data and re-clustered. Cells from each sample (Supplementary Figure 2A) and model (Supplementary Figure 2B) were visualized to confirm the batch correction. Although some cells showed biased distribution in a few clusters of the integrated data, suggesting different cell numbers in some samples in these clusters, most of the cells were evenly distributed (Supplementary Figures 2C,D).

We next determined the transcript of basal marker Krt5 and the luminal marker Krt18. Notably, all the cancer cells from the four models showed high Krt18 expression, and only some cells from the BRCA1 model, a basal-like breast cancer model, showed a detectable Krt5 transcript (Supplementary Figure 2E). We also checked the transcriptional levels of ER genes, PR genes, and HER2 (Erbb2) and found rare transcripts in the cancer cells of all the tumor models, even in the NEU model (Supplementary Figure 2F), suggesting that these were good models to mimic TNBC.

As the control of the tumor cells, we then accessed another scRNA-seq dataset on normal mouse mammary epithelial cells at four developmental stages: nulliparous, mid-gestation, lactation, and post involution (Bach et al., 2017). Similarly, sample integration and QC were performed on these data (Supplementary Figures 3A,B). After clustering, cells from different samples and developmental stages, as well as the cell clusters, were visualized by UMAP (Supplementary Figures 3C–E). The major cell types, including epithelial cells, mesenchymal cells, immune cells, and endothelial cells, were identified by checking canonical cell-type marker genes (Supplementary Figures 3F–H). The epithelial cells were then subset and re-clustered, and epithelial cells from different stages (Supplementary Figures 4A,B) showed biased distribution in the clusters (Supplementary Figures 4C,D), suggesting that the gene profiles in different stages differed significantly, confirming the findings in the original study (Bach et al., 2017). The basal marker Krt5 and the luminal marker Krt18 were used to confirm this (Supplementary Figures 4E,F). It was clear that the epithelial cells in nulliparous and mid-gestation stages were composed of both Krt18⁺ and Krt5⁺ cells, and in the lactation stage, the epithelial cells were mostly Krt5⁺ cells, while in the post involution stage, the epithelial cells were mainly Krt18⁺ cells (Supplementary Figures 4E,F).

The cancer cells were subset from the tumor models, and the normal epithelial cells were subset from normal mammary datasets, and then these cells were merged (Figure 1A). Due to the gene profile differences, the cancer cells and normal epithelial cells showed distinct clusters (Figure 1B). The integrated data were clustered, and 17 major cell clusters were identified (Figure 1C). Most of the clusters showed specific expression of either the luminal markers Krt18 and Krt8 or basal markers Krt14, Krt5, and Acta2 (Figure 1D).

To determine the gene profiles of the cells of different developmental stages/models, we performed comparative analysis on the scRNA-seq data, and the top 100 genes were presented by heatmap, suggesting distinct gene profiling patterns in each developmental stage/model (Figure 1E). Further comparative analysis was used to identify the profiling difference of normal epithelial cells and cancer cells. The top 500 genes in the total normal epithelial cells and total cancer cells were visualized by heatmap (Figure 1F). The top 20 upregulated genes in normal epithelial cells (Figure 1G) and in cancer cells (Figure 1H) were visualized by violin plots. Most of the top genes identified in normal epithelial cells, including Acta2 (Weymouth et al., 2012), Mylk, Cnn1 (Bresson et al., 2018), and multiple Krt genes (Kendrick et al., 2008), were associated with normal development and functions of mammary gland, such as maintenance of lactation, mammary stem cell function, and lineage commitment determination. Notably, two function-unidentified ribosome protein genes, Rpl13a and Rpl29, were significantly downregulated in cancer cells (Figure 1G and Supplementary Figure 5A). Of the top genes upregulated in the cancer cells, most were reported to be breast cancer-promoting genes. Specifically, AGTG1 is one of the recurrent mutated genes in breast cancer (Ciriello et al., 2015). Lalba is reported to not be expressed in healthy tissues except during lactation, while it is expressed at high levels in a subset of human breast cancers, especially TNBC (Tuohy, 2014). Higher SPP1 gene expression in primary tumors was found to be associated with risk of recurrence in ER⁺ breast cancer among patients with endocrine treatment (Gothlin Eremo et al., 2020). Modification of Sumo2 was associated with the progression and metastasis of breast cancer (Subramonian et al., 2014). Other genes were associated with cell adhesion and migration of breast cancer cell [DUOXA1 (Ostrakhovitch and Li, 2010) and PGAM1 (Zhang et al., 2017)], prognosis of breast cancer [PGK1 (Fu et al., 2018) and KLK10 (Kioulafa et al., 2009)], or clinical outcome in breast cancer [RBM3 (Jogi et al., 2009) and AQP5 (Lee et al., 2014)]. Another gene was Rpl15, which was recently reported to markedly increase metastatic burden in TNBC mouse models (Ebright et al., 2020). This was further confirmed in our study by the data that an increased expression of Rpl15 gene was found in the cancer cells. A more interesting gene was Rpl27a, which also showed significantly increased transcript in cancer cell (Figure 1H and Supplementary Figure 5A) and whose functions had rarely been studied in diseases or cancers.

To better clarify the expression level of Rpl27a, the transcripts were visualized by both violin plots and UMAPs (Figure 2A and Supplementary Figure 5A). The transcripts of Rpl27a in split violin plots (Figure 2B) and UMAPs (Supplementary Figure 6) further confirmed the significantly upregulated transcription in cancer cells of all the cancer models compared to normal epithelial cells in all the mammary developmental stages. Consistently, Rpl15 was also significantly upregulated and Rpl29 and Rpl13a were downregulated (Figures 2A,B and Supplementary Figures 5A, 6). To further quantify the expression of Rpl27a, we calculated the percentage of the cells with Rpl27a expression at Rpl27a > 3. Notably, percentages of cells with Rpl27a expression among the four tumor models were comparable but were significantly higher than those of normal mammary epithelial cells of developmental stages (Figure 2C). This was also true for Rpl15 (Figure 2D). On the contrary, percentages of Rpl29-expressing and Rpl13a-expressing cells were notably decreased in cancer models (Supplementary Figures 5B,C). These data confirmed that Rpl27a was potentially a diagnosis biomarker for TNBC in the mouse models, but its expression levels in human breast cancer and its functions and underlying mechanisms needed to be further determined.

To validate the expression of Rpl27a in human TNBC tumor cells, we first accessed a published scRNA-seq data on primary human TNBC tumors (Karaayvaz et al., 2018). Similar to the analysis of mouse scRNA-seq data, we performed QC to remove the low-quality cells (Supplementary Figures 7A,B). High-quality cells were then clustered (Supplementary Figure 7C), and major cell types were identified by canonical cell-type marker expression (Supplementary Figures 7D,E). The cancer cells were subset for further analysis (Supplementary Figure 7F). We also extracted the scRNA-seq data on human metastatic TNBC cells to the brain from another dataset. Two batches of data were integrated, and QC was performed to remove the low-quality cells (Supplementary Figures 8A,B). Retained cells were clustered (Supplementary Figure 8C), and major cell types were defined (Supplementary Figures 8D,E). Cancer cells were then subset for further analysis (Supplementary Figure 8F).

The purified primary and metastatic cancer cells were integrated and clustered (Figure 3A), and three major cell clusters were identified (Figure 3B). Quantitatively, primary cancer cells shared equal percentage in cluster 0 and cluster 2 and had lower percentage in cluster 1 (Figure 3C). To determine the lineage differentiation potentials of the cancer clusters, we loaded the cells into Slingshot packages, packages that provided unified interface to dozens of different trajectory inference methods via docker containers, for pseudotime analysis. The cancer cell clusters were matched to the clusters defined above, and the development went along a cluster 1–cluster 0–cluster 2 axis (Figure 3D). Analysis on inferred trajectory and principal curves demonstrated smoothed representations of each subcluster (Figure 3E). Slingshot analysis revealed a pseudotime starting from cluster 1 and ending at cluster 2 (Figure 3F), which is consistent with the quantification analysis that primary cells contributed higher-percentage cells to cluster 1 and lower-percentage cells to cluster 2 compared to metastatic cells (Figure 3C, bottom panel). This analysis suggested that most metastatic cells (cluster 2) were mostly derived from less metastatic cells (clusters 0 and 1).

Comparative analysis revealed gene expression patterns of the three clusters, which was presented by the top 100 genes of each cluster (Figure 3G). Gene profiles visualized by the top 500 genes confirmed the specific gene expression patterns of both primary and metastatic cancer cells (Figure 3H). We then visualized the top 20 upregulated genes of both primary and metastatic cancer cells and showed their expression by violin plots (Figures 3I,J). Many of the genes upregulated in primary cancer cells showed some potential association in breast cancer progression and metastasis inhibition. In both genetically engineered mouse models and xenograft models, lncRNA MALAT1 overexpression inhibits, while MALAT1 deficiency induces breast cancer metastasis, which is reversed by re-expression of MALAT1, suggesting that MALAT1 is a metastasis-suppressing lncRNA (Kim et al., 2018). SRSF3 and HNRNPH1 regulate a splicing hotspot of HER2 in breast cancer cells (Gautrey et al., 2015). Another lncRNA KCNQ1OT1 is correlated with human breast cancer cell development through inverse regulation of miR-145 (Feng et al., 2018). TRA2A promotes paclitaxel resistance and tumor progression in TNBC via regulating alternative splicing (Liu et al., 2017). Notably, majority of the genes upregulated in metastatic cancer cells were ribosome protein genes (Figure 3J), suggesting that these genes were closely related with TNBC cell metastasis.

Among the ribosome protein genes, the RPL27A and RPL15 genes were upregulated in metastatic cancer cells with significant differences, which was consistent with the data in mouse breast cancer models (Figures 4A,C). The split UMAP visualization of RPL27A and RPL15 expressions in primary and metastatic cancer cells confirmed these differences (Figures 4B,D). As RPL15 was a reported pro-metastasis gene in breast cancer (Ebright et al., 2020), we next checked the transcript colocalization of RPL27A with RPL15 in primary, metastatic, and integrated data. Surprisingly, RPL27A was well colocalized with RPL15 only in metastatic cancer cells; however, in primary cancer cells, the transcripts of the two genes were low and rarely colocalized (Figure 4E). Although Rpl29 and Rpl13a were downregulated in cancer cells of mouse models compared to normal mammary epithelial cells, these two genes were significantly upregulated in metastatic cancer cells compared to primary cancer cells (Supplementary Figure 9). These inconsistent findings suggested diverse functions of these genes or an inconsistent phenotype between mouse models and human tumor data (Supplementary Figure 9). Notably, the transcript levels of these four genes in total cells were also upregulated in metastatic TNBC tissues (Supplementary Figure 10A).

To further confirm the potentials of RPL27A in regulating cancer cell metastasis, an in vitro model was generated on a human TNBC cell line, MDA-MB-231. RPL27A knockdown was performed on MDA-MB-231 cells by shRPL27A adenovirus transfection, and the knockdown efficiency was confirmed by quantitative real-time PCR (qRT-PCR) and western blot (Figures 5A–C). The cells were then used for migration and invasion assays, and significantly reduced cell numbers in both migration and invasion assays were confirmed by cell quantification (Figures 5D–G). These data further suggested that RPL27A was required for TNBC cell metastasis. More interestingly, RPL27A knockdown also caused significant downregulation of many other ribosomal genes (Supplementary Figure 10B).

As many of the ribosomal protein genes were in the top gene list of the metastatic human TNBC cells, we proposed that the ribosomal protein genes were deregulated in TNBC metastasis. We then listed all the significantly differentially expressed ribosome protein genes and visualized them by heatmap (Figure 6A), and among them RPL27A, RPL15, RPL13A, and RPL29 were the top genes. We also separated the RPL genes and RPS genes and scored their expression by UMAPs (Figures 6B,C). Both the scores of RPL genes and RPS genes were significantly higher in metastatic cancer cell compared to those of primary cancer cells (Figures 6B,C). We next input these genes and their corresponding expression fold changes and p-values into IPA, and the top ingenuity canonical pathways and upstream regulators were visualized by bar plots (Supplementary Figures 10C,D). The graphic regulatory network based on the IPA revealed that EIF2 signaling was the core regulating signaling pathway (Supplementary Figure 10E).

Epithelial–mesenchymal transition (EMT) is one of the major characteristics of breast cancer metastasis. EMT programs were believed to reflect a loss of epithelial gene expression signatures and morphologies that gave way to those associated with mesenchymal cells and their enhanced migratory and invasive behaviors (Gooding and Schiemann, 2020). To check the EMT activities of the primary and metastatic cancer cells, the transcriptional levels of canonical EMT-related signature genes, including epithelial feature genes (CDH1, KRT8, KRT18, CTNNB1, and ZEB1), mesenchymal feature genes (VIM, S100A4, FN1, and LAMA5), and EMT signature genes (SNAI1, SNAI2, ITGB6, TGFB1, and TGFB2) were visualized and compared between primary and metastatic cancer cells (Figure 7A). Significantly decreased epithelial signatures and increased mesenchymal and EMT signatures were found in metastatic cancer cells. Scoring of the epithelial signatures and mesenchymal signatures on primary and metastatic cancer cells further confirmed decreased epithelial scores and increased mesenchymal scores in metastatic cancer cells (Figures 7B,C).

To determine the major cancer metastasis-related signaling pathways, the differentially expressed genes between primary and metastatic cancer cells were input into IPA. The top 1 ingenuity canonical pathway was EIF2 signaling, and this (Figure 7D), together with the IPA on ribosome protein genes, further confirmed that EIF2 signaling was the key driver of the TNBC metastasis in human. To further confirm that, we checked the expression of the major component genes of EIF2 signaling and found upregulation of multiple genes, including EIF2AK1 (HRI protein coding gene), EIF2AK2 (PKR protein coding gene), EIF2AK4 (GCN2 protein coding gene), EIF2B5, and AKT2 (Figure 7E; Rios-Fuller et al., 2020). Other ingenuity canonical pathways included regulation of eIF4 and p70S6K signaling and mTOR signaling (Figure 7D), both of which were closely associated with TNBC biological activities (Madden et al., 2014; Wang et al., 2017; Zagorac et al., 2018). Graphic regulatory networks integrating ingenuity canonical pathways and upstream regulators of the metastatic cells gave a core pathway of cell viability of breast cancer cells, which was closely related with EIF2 signaling and regulators of RB1, TP53, and MYC (Figure 7F).

The salubrinal derivative Sal003 is the most commonly used eIF2α inhibitor II (Robert et al., 2006; Baltzis et al., 2007). To better validate the functions of EIF2 in TNBC metastasis, we inhibited EIF2 signaling in the MDA-MB-231 cells with Sal003 at different doses. Consistently, both the mRNA levels of RPL27A and RPL15 were significantly decreased after Sal003 treatment at a dose-dependent manner (Figures 8A,B), and this is also true for many other ribosome-related genes we detected (Supplementary Figure 10F). More importantly, Sal003 treatment dramatically attenuated MDA-MB-231 cell migration and invasion (Figures 8C–F), suggesting an attractively reduced metastasis potential in cancer cells when blocking EIF2 signaling. All these findings supplied basic foundations for further studies and potential therapeutic target development of TNBC and were likely to attract a broad audience.