TILRR (FREM1 isoform 2) is a prognostic biomarker correlated with immune infiltration in breast cancer

FREM1 mRNA expression levels in human breast cancer tumors

To investigate the role of FREM1 gene expression in cancer, we determined the FREM1 mRNA levels in tumor and normal adjacent tissues of 20 cancer types. Data obtained from Oncomine was analyzed using the following threshold criteria: 2-fold change, P value < 0.0001 and a gene rank of 10%. Downregulation of FREM1 gene expression was found in breast cancer, ovarian cancer and pancreatic cancer tissues. In breast cancer tissues, 9 out of 43 samples met the threshold criteria in 4 out of 10 datasets. In the other two cancer types, only a single event reached threshold criteria (Figure 1A).

Further evaluation of FREM1 expression in multiple malignancies was carried out using the gene expression data from The Cancer Genome Atlas (TCGA). Differential FREM1 expression between the tumor and adjacent normal tissues is illustrated in Figure 1B. In most types of cancer, FREM1 expression is significantly lower (P value < 0.001) than in adjacent normal tissues (Figure 1B).

FREM1 expression in breast cancer was further investigated using TCGA and Gene Expression Omnibus (GEO) data (Figure 1C and 1D). FREM1 expression in tumor tissues was 1.25-fold (GEO) and 4.18-fold (TCGA) lower than in normal adjacent tissues. These findings were confirmed in paired patient biopsy samples via Q-PCR (N = 36, P = 0.0336) (Figure 1D and 1E). In tumor biopsy samples, FREM1 expression was reduced by as much as 6.16-fold (P < 0.001). In 17 samples (47.2%), expression was decreased more than 2-fold (Figure 1E). We expanded our analysis of FREM1 expression to various subtypes of breast cancer (Figure 1G and 1H). This analysis compared normal tissue (n = 120/n = 75) with tissues from the following breast cancer subtypes: Basal-like (n = 168/n = 75), HER2-enriched (n = 78/n=24), Luminal-like A (n = 493/n = 89) and Luminal-like B (n = 194/n = 49). These data showed decreased FREM1 expression in Basal-like (6.07-fold/1.11-fold with mean), HER2-enriched [6.79-fold/1.09-fold with mean (P = 0.0279)], Luminal-like A (3.46-fold/1.08-fold with mean) and Luminal-like B (7.54-fold/1.11-fold with mean) tissues. Hence, FREM1 expression is downregulated in breast cancer tissue.

FREM1 transcription level is correlated with survival and progression in breast cancer patients

We further analyzed TGCA data to determine if there is a relationship between FREM1 expression and overall survival (OS) or disease-free survival (DFS) in breast cancer patients. Patients were divided into high- or low-level groups over the median value of FREM1 expression in breast cancer tissues. Patients with high FREM1 expression in their tumors experienced a prolonged OS and DFS compared with those who with low FREM1 expression (Figures 2A and 2B). Additionally, KM analysis revealed that HER2-positive BRCA patients in the low-level FREM1 expression group generally demonstrated shorter DFS (Figures 2C and 2D). To confirm the clinical significance of downregulated FREM1 expression in breast cancer patients, the correlation between survival rate and FREM1 transcription was determined. Survival analysis of GEO data indicated that low-level FREM1 expression was correlated with reduced OS and DFS (Figure 2E and 2F), and the same scenario was shown in HER2-positive BRCA patients (Figure 2G and 2H). This data was consistent with the results obtained from TCGA cohort. These results suggest that reduced FREM1 transcription influences breast cancer tumor progression and is associated with shortened patient survival.

Downregulation of the FREM1 protein in primary human breast cancer tissues correlates with disease progression

To determine the clinical significance of the FREM1 protein, we performed representative immunohistochemical (IHC) staining in primary human breast cancer tumor and adjacent normal tissue (N = 47). FREM1 staining in human primary breast cancer was scored as follows: +++, high; ++, moderate; +, weak; −, negative. FREM1 expression was significantly decreased in tumors, as indicated by lower IHC staining scores (P < 0.001) (Figure 3A and 3B). Among the 47 breast cancer samples, FREM1 staining was high in 4 samples (8.5%), moderate in 22 samples (46.8%) and weak or undetectable in 21 samples (44.7%) (Figure 3C and 3D). The AOD (average optical density) quantified data set was plotted in Supplementary Figure 1.

From the clinical information provided with the tissue chip, the survival rate was analyzed. High FREM1 expression correlated with a prolonged OS and DFS compared with low FREM1 expression, but the statistical significance of the DFS prolongation was weak (P = 0.0509) (Figure 3E and 3F). When data was re-analyzed only using HER2-enriched samples, the correlation between FREM1 expression and survival rate was robustly significant in both OS (P = 0.0006) and DFS (P = 0.0035) (Figure 3G and 3H). Thus, FREM1 downregulation in breast cancer tissue, especially in HER2-enriched samples, was correlated with shorter survival times. This finding from IHC staining is consistent with our database analysis.

TILRR is the clinically relevant isoform of FREM1 in breast cancer

Previously, TILRR has been identified as isoform 2 of the FREM1 gene [3]. We wanted to determine if TILRR is the dominate isoform that is downregulated in breast cancer. In the Ensembl and Gene Expression Profiling Interactive Analysis 2 (GEPIA2) databases, three different FREM1 isoforms are described. In these isoforms, ENST00000422223.6 and ENST00000380880.3 encode a protein with 2179 amino acids, and the much shorter version of ENST00000380894.5 encodes a protein with 715 amino acids. ENST00000380894.5 is the FREM1 isoform 2 transcript, which we referred to previously as TILRR (Figure 4A). The correlation between breast cancer survival rate and FREM1 usage was investigated using the GEPIA2 database. In Figure 4B, the violin-plot and bar-plot panels present the expression distribution and usage of each FREM1 isoform. TILRR, or FREM1 isoform 2, shows convincing cancer specificity with strong usage and a significantly low hazard ratio (HR).

Using a TILRR-specific primer (Figure 4C), TILRR expression in 30 paired breast cancer samples was investigated using Q-PCR. The data showed a 4.39-fold decrease (P = 0.0274) in expression in tumor tissue compared with normal tissue (Figure 4D and 4E), which confirmed that TILRR is downregulated in breast cancer tissues. The protein expression identity of FREM1 and TILRR was verified using four paired breast cancer tissues (sample numbers: 160, 172, 109 and 164) (Figure 4F and 4G). The FREM1 commercial polyclonal antibody blotting pattern was identical to that of the TILRR polyclonal antibody, which detected a 75 KDa peptide with an approximately 4.76- and 9.09-fold reduction in expression in tumor tissues. A similar finding was discovered using TILRR monoclonal probes (data not shown). These results confirm that downregulation of TILRR is associated with breast cancer prognosis.

Analysis of TILRR transcript level associated with tumor-infiltrating lymphocytes

As shown above, TILRR is the dominant isoform of FREM1 expressed in breast cancer tissues. A prior study found that TILRR is related to monocyte infiltration in atherosclerosis plaque development [4]. Moreover, the lymphocyte-specific immune (LYM) recruitment metagene signature is related to tumor infiltration by lymphocytes, and it is associated with a favorable prognosis in breast cancer. To investigate the role of TILRR in the mediation of immune cell infiltration, the correlation between TILRR expression and LYM was evaluated in TCGA breast cancer data. The LYM metagene sets PTPRC (CD45), CD53, LCP2 (SLP-76), LAPTM5, DOCK2, IL10RA, CYBB,CD48, ITGB2 (LFA-1) and EVI2B were positively correlated with expression of FREM1 (Spearman r = 0.42 P = 3.5e-47; Pearson r = 0.3 P = 0) (Figure 5A).

The estimation of stromal and immune cells in malignant tumor tissues using expression data(ESTIMATE) algorithm can calculate the gene expression signature in tumor cells and normal, tumor-associated epithelial and stromal cells, immune cells and vascular cells [5]. Stromal cells are thought to have essential roles in tumor growth, disease progression and drug resistance, and ESTIMATE immune scores could serve as an indicator for immunotherapy response [5]. Tumor purity and the expression of TILRR were analyzed using the ESTIMATE algorithm to investigate in GEO (N = 237) and TCGA (N = 1097) data. The results showed that the expression of TILRR was negatively correlated with tumor purity, but positively correlated with immune score and stromal score (Figure 5B). This correlation was confirmed with the 50 highest and 50 lowest TILRR-expressing tumor tissue samples (Figure 5C and 5D).

To understand the underlying mechanism of FREM1 in breast cancer, we analyzed differentially expressed data from the 50 highest and 50 lowest FREM1-expressing tumor tissue samples. The gene expression enrichment of several signal pathways was analyzed using gene set enrichment analysis (GSEA). TILRR is positively correlated with Toll-like receptor signaling and cytokine-to-cytokine receptor interactions (Figure 5E), which is consistent with our previous finding [3]. We also found that leukocyte transendothelial migration, chemokine signaling, T-cell receptor signaling and JAK-STAT signaling are associated with the migration and infiltration of immune cells (Figure 5E). To verify these pathways, RNA-seq analysis of the BT474 cancer cell line overexpressing TILRR was performed using an empty GFP construct as a control. The analysis showed that TILRR overexpression not only potentiates IL-6 and CXCL8 immune- and inflammation-induced gene expression, but it also induces the expression of other cytokine genes such as CCL5/2 and chemokine genes such as CXCL10 (Figure 5F). JAK-STAT signaling pathways and interferon-stimulated genes (further detailed analysis will be published elsewhere) were also upregulated (Figure 5F). This upregulation of pro-inflammatory secretory factors and pathway activation indicated that enhanced TILRR expression might cause long distal immune recruitment. Q-PCR was performed to confirm RNA-seq detection among expression of CXCL8, CXCL10, CXCL11,MX1 and the interferon stimulation gene, ISG15 (Figure 5G).

Furthermore, we analyzed the 50 highest TILRR-expressing and 50 lowest TILRR-expressing tumor tissue samples (TCGA database) to determine the relationship between the expression of CXCL10 and CXCL11 and the expression of TILRR in cancer tissues. The results showed that high expression of TILRR positively correlates with CXCL10 and CXCL11 expression (Supplementary Figure 2). These findings strengthened the hypothesis that TILRR expression is associated with breast cancer progression and prognosis, likely through signaling pathways that regulate the distal recruitment of immune cell infiltration.

Expression of TILRR is correlated with immune cell infiltration in breast cancer

Increasing evidence suggests that tumor-infiltrating immune cells can be an indicator in the clinical analysis of tumor samples [6, 7]. Gene expression profiling of heterogeneous cell populations in cancer tissue, including tumor-infiltrating lymphocytes, serves as an independent predictor of survival in prognostic cancer models [8, 9]. Thus, we used the TIMER database to evaluate the correlation between FREM1 mRNA expression and six different infiltrating immune cell types (B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils and dendritic cells) in different subtypes of breast cancer. Data showed that the TILRR transcription level was inversely related to the purity of tumor tissue in breast cancer subtypes (Figure 6A–6D, first line on the left). The low heterogeneity in tumor tissue correlated with high expression of FREM1 in all but one subtype, which is consistent with the positive association between TILRR expression and level of immune cell infiltration (Figure 6). The exception was macrophage cells of the HER2 subtype (r = -0.107, P = 4.23e-01). The infiltration levels of CD8+ and CD4+ cells were significantly positively correlated with TILRR transcription. The level of CD8+ T-cell infiltration with TILRR expression in BRCA (r = 0.379, P = 16e-34), BRCA-Basel (r = 0.275, P = 2.10e-03), and BRCA-HER2 (r = 0.406, P = 1.2e-03), and BRCA-Luminal (r = 0.48, P = 7.89e-28) were all significantly strong, respectively, which is similar to the CD4+ T-cell infiltration levels (BRCA: r = 0.369, P = 2.07e-32; Basal: r = 0.181, P = 4.57e-02; HER2: r = 0.595, P = 8.37e-07; luminal: r = 0.42, P = 2.44-24e). On the other hand, in the tumor tissue, the infiltration level from B cells, macrophages, neutrophils and dendritic cells correlated more weakly with TILRR transcription in subtypes of breast cancer (Figure 6). To verify the expression of TILRR and the infiltration of immune cells, the 50 highest TILRR-expressing and 50 lowest TILRR-expressing tumor tissue samples were selected (TCGA database) to analyze the expression of marker genes for CD8+, CD4+ and T cells (general). The results showed that CD8A, CD8B was expressed in tissues with high TILRR expression. (Figure 6E). Analysis of GEO data generated similar results (Supplementary Figure 3). To further verify this result, we selected the 5 samples with the highest TILRR expression and the 5 with the lowest TILRR expression patients shown in Figure 4D to evaluate the expression of marker genes in CD4- and CD8-expressing cells. Among these 10 samples, the expression of CD8A, CD8B and CD4 in TILRR high- and low-expressing samples were positively correlated (Figure 6F). We also investigated the correlation between TILRR expression and biomarkers from other immune cells. Strong associations were found between NK cells, Treg cells, T-cell exhaustion and TILRR expression (Supplementary Table 1). Consequently, it is likely that TILRR expression in breast cancer tissue is involved in immune cell recruitment.

FREM1 gene expression of BRCA in the TCGA and GEO databases

The expression level of the FREM1 gene in various types of cancers was identified in the Oncomine database (https://www.oncomine.org/resource/login.html) [29]. The threshold was determined according to the following values: P value of 0.001, fold change of 2. Differential expression module of TIMER (https://cistrome.shinyapps.io/timer/) was used to analyze the FREM1 differential expression between tumor and adjacent normal tissues of various TCGA tumors [30, 31]. TCGA breast cancer data was downloaded from the Xena browser (https://xenabrowser.net/datapages/) [32]. GSE21653 data was downloaded from the GEO database, which contained 266 early cancer patients [33, 34].

Immune cell infiltration analysis

The correlation of FREM1 expression level with a LYM metagene signature [PTPRC (CD45), CD53, LCP2 (SLP76), LAPTM5, DOCK2, IL10RA, CYBB, CD48, ITGB2 (LFA-1) and EVI2B] was explored via Gene Expression Profile Interactive Analysis (GEPIA2; http://gepia2.cancer-pku.cn/#index) [35–37]. The stromal score, immune score and tumor purity of breast cancer (BRCA) tumor samples from the GSE21653 database (N = 237) and the TCGA database (N = 1097) were calculated by using the R 3.6.0 ESTIMATE package (1.0.13) [38]. The gene module of TIMER was used to evaluate the correlation of FERM1 expression with immune cell infiltration (including B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages and dendritic cells) in BRCA, which determined the purity-corrected partial Spearman’s correlation and statistical significance. The correlation module of TIMER was used to calculate the Spearman’s correlation and statistical significance between FREM1 and immune cell marker genes. The top 50 highest and lowest FREM1-expressing TCGA tumor samples were used to conduct GSEA and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotations analysis [39].

RNA isolation and real-time PCR

Total mRNA of breast cancer tissue (SYSUCC, Sun Yat-sen University Cancer Center, Guangzhou, China) was extracted using TRIzol reagent (Invitrogen, #1556018) according to the manufacturer’s protocols. The quantitative, real-time PCR using SsoAdvanced^TM Universal SYBR Green Supermix (Bio-Rad, #1725274) was performed in an ABI StepOnePlus (ABI, #1725274). FREM1 primer sequences: forward primer, 5’-AGAGCCCTGCCTGTGGTAAC-3’; reverse primer, 5’-GAAGGGGAATGCAAGAGTGTGATA-3’. TILRR-specific primer 5’-GCCTTGCCTCTCTTTACCAGAT-3’; reverse primer, 5’-GAGTGCCGATAGGCCACAT-3’. Relative gene expression was normalized to glyceraldehyde 3-phosphate dehydrogenase (GAPDH) (forward primer, 5’-GGAGCGAGATCCCTCCAAAAT-3’; reverse primer, 5’-GGCTGTTGTCATACTTCTCATGG-3’) expression and was analyzed using the 2^−ΔΔCT method.

Immunohistochemistry (IHC)

Immunohistochemistry was performed on tissue section microarrays (Zhuolibiotech, #ZL-Brc Sur122). The staining procedure included heat-induced epitope retrieval using 0.1 M sodium citrate (pH 9.0, heated by microwave to 90-95 °C, 3 times for 5 minutes each), incubation with primary antibody at 4 °C overnight. Signal detection was performed using an IHC detection kit (Gene Tech, #GK500710). Microscopy of the immunostaining included an initial pre-screen at low power (4X) to identify regions with a technically optimal staining result. Subsequently, detailed analysis at high-power (40X) was performed to evaluate the staining according to routine algorithms employed in tumor diagnostics. Scoring of FREM1 staining were evaluated using the semi-quantitative immunostaining score (ISS) method by pathologist. The immunostaining score was defined as 0 – 3 (range: +++/3, high; ++/2, moderate; +/1, weak; 0, negative.) [40]. The median score was used as cut off for classification of patients into high- and low- expression groups. Semi-quantitative analysis of the IHC images was conducted by Image-J, by which integral optical density (IOD) and the area were collected. Then, average optical density (AOD) was calculated as IOD/area, which represented the staining intensity [41].

Western blot

Total proteins were extracted from breast cancer and adjacent cancer tissue (SYSUCC, Sun Yat-sen University Cancer Center, Guangzhou, China) by using RIPA buffer (Beyotime, #P0013B) and passed the tissue through a 21-gauge needle more than 30 times. The cell lysate was centrifuged at 4 °C and 17000 g for 15 minutes to eliminate large aggregates. A BCA Protein Assay Kit (Tiangen, #PA115) was used to evaluate the cell lysate protein concentrations. Approximately 25-50 μg total protein, denatured with the standard SDS-sample buffer per lane, was loaded for electrophoresis on 10% pre-cast mini-polyacrylamide gels (GenScript SurePAGE, Bis-Tris, 10 cm x 8 cm gels). The gel containing proteins was then transferred to a PVDF membrane (BioRad). The membrane was blocked with 5% milk in 1×TBST for 2 hours, and incubated with the FREM1 primary antibody (1:1000 rabbit polyclonal, Proteintech, #13086-1-AP) and TILRR (1:2000 rabbit polyclonal custom, Genesript) at 4 °C overnight. Before incubating the membrane with a secondary antibody (1:10000 dilution, Jackson ImmunoResearch, #111-035-003) for 1 hour, it was washed the with 1×PBST 5 times at room temperature. The blots were detected by chemiluminescence (Bio-Rad). The band intensity was quantified by ImageJ (https://imagej.nih.gov/ij/download.html).

Cell culture

Cells were cultured (BT474, purchased from Procell Life Science and Technology Co., Ltd. CL-0040) at the density of 1×10⁶/well in RPMI 1640 Medium (Gibco, #C11875500BT) with 20% FBS and insulin. Cells were incubated at 37 °C and 5% CO₂. The day after seeding, cells were transfected with HA-GFP or HA-TILRR-T2A-GFP plasmid using the Lipofectamine 3000 Transfection Reagent (Invitrogen, #L3000-015). Cells were incubated at 37 °C and 5% CO₂ for 24 hours, at which time culture medium was replaced. Forty-eight hours after transfection, cells were detached by 0.25% trypsin for 3 minutes, and then collect the cell for cell sorting (Flow Cytometer, Aria II, BD).

RNA-seq

Trizol lysate was used to lyse selective cells and extract RNA from the cells. To generate sequencing libraries, 1 ug RNA per sample was used (NEBNext UltraTM RNA Library Prep Kit, Illumina; NE, USA). Library quality was assessed on the Agilent Bioanalyzer 2100 system. Sequencing was performed on the Illumina Novaseq platform, and 150 bp paired-end reads were generated. Raw reads were aligned to a reference genome (UCSC GRCh38/hg38) using Bowtie 2 (2.2.5). Gene expression was quantified using RSEM v1.1.22.

Statistical analysis

OS and DFS curves were generated by Kaplan-Meier survival analysis using SPSS 17.0 software. The results generated in Oncomine are displayed with P values, fold changes, and ranks. The results of the Kaplan-Meier plots, PrognoScan, and GEPIA2 are displayed with HR and P or Cox P values from a log-rank test. The correlation of gene expression was evaluated by Spearman’s correlation and statistical significance, and the strength of the correlation was determined using the following guide for the absolute value: 0.00–0.19, very weak; 0.20–0.39, weak; 0.40–0.59, moderate; 0.60–0.79, strong; 0.80–1.0, very strong. P values < 0.05 were considered as statistically significant.