Prognostic implications of the expression levels of different immunoglobulin heavy chain-encoding RNAs in early breast cancer

Abstract

The extent and composition of the immune response in a breast cancer is one important prognostic factor for the disease. The aim of the current work was to refine the analysis of the humoral component of an immune response in breast tumors by quantifying mRNA expression of different immunoglobulin classes and study their association with prognosis. We used RNA-Seq data from two local population-based breast cancer cohorts to determine the expression of IGJ and immunoglobulin heavy (IGH) chain-encoding RNAs. The association with prognosis was investigated and public data sets were used to corroborate the findings. Except for IGHE and IGHD, mRNAs encoding heavy chains were generally detected at substantial levels and correlated with other immune-related genes. High IGHG1 mRNA was associated with factors related to poor prognosis such as estrogen receptor negativity, HER2 amplification, and high grade, whereas high IGHA2 mRNA levels were primarily associated with lower age at diagnosis. High IGHA2 and IGJ mRNA levels were associated with a more favorable prognosis both in univariable and multivariable Cox models. When adjusting for other prognostic factors, high IGHG1 mRNA levels were positively associated with improved prognosis. To our knowledge, these results are the first to demonstrate that expression of individual Ig class types has prognostic implications in breast cancer.

Keywords: Breast cancer; Prognostic markers; Tumour immunology.

Conflict of interest statement

Competing interestsThe authors declare that they have no competing interests.

© The Author(s) 2020.

Figures

Fig. 1. Expression and correlation of RNAs encoding IGJ and Ig heavy chains of different types.

a, b Box plot of the log2 expression values of RNA encoding IGJ and the constant region of the heavy chain of different Ig classes. The data are from SCAN-B (a) and Cohort270 (b). The center line marks the median; box limits mark the upper and lower quartiles; whiskers mark 1.5× interquartile range; points mark outliers beyond this mark. Correlation matrix of the log2 expression of the indicated mRNAs in SCAN-B (c) and Cohort270 (d). The color indicates the level of the Pearson’s correlation coefficient. e, f The microenvironment cell population counter was used to obtain quantitative data for stromal cell types. Pearson’s correlation coefficient between the cell type values and log2 IG mRNA levels were calculated for SCAN-B (e) and cohort270 (f). The color indicates the levels of the correlation coefficient. g The Spearman’s correlation between indicated IG mRNAs and the amount of lymphocytes in the tumor, estimated on tumor section.

Fig. 2. Enrichment of gene ontology sets among genes correlating with IGHA2 and IGHG1 mRNAs.

All mRNAs in the SCAN-B cohort were analyzed for correlation with IGHA2 and IGHG1 mRNA expression. The top 100 correlating genes for each IG mRNA were selected for enrichment analysis. The enrichment was analyzed performed using the Biological Process sets defined by the Gene Ontology Consortium. Fisher’s test was used to evaluate enrichment. The graph shows the −log10 of the p value from the Fisher’s tests of indicated Gene Ontology Biological Process gene sets for the top 100 IGHA2- (blue) and IGHG1- (red) correlating genes. Results are shown for gene sets with a p value < 10–15 for either IGHA2- or IGHG1-correlating genes.

Fig. 3. Expression of Ig heavy chain-encoding RNAs in relation to clinico-pathological parameters.

Box and scatter plots of the log2 expression levels of indicated RNA species in relation to ER status, HER2 amplification, NHG grade, PAM50 subtypes, and age at diagnosis for the SCAN-B cohort. The p values were calculated with the t-test comparing positive with negative status (ER and HER2) and grade 1 versus grade 3. For age the p value was extracted from linear regression modeling using tlm function in R. The PAM50 subtypes are basal (B), HER2-enriched (H2), luminal A (LA), luminal B (LB), and normal-like (N). In box plots the center line marks the median; box limits mark the upper and lower quartiles; whiskers mark 1.5x interquartile range; points mark outliers beyond this mark.

Fig. 4. Association of IG mRNA expression with overall survival.

Breast cancers from four data sets were grouped based on if the expression of IGHM, IGHA2, IGJ, and IGHG1 mRNA was larger (red curves) or smaller (black curves) than the median expression. The figures display Kaplan–Meier curves for the indicated mRNAs (top titles) and cohorts (titles to the left) using overall survival as end point. The follow-up time on the x-axis is indicated in years. The p values were estimated with the log rank test. The y-axes indicate years after diagnosis and number at risk in the groups.

Fig. 5. Association of IG mRNA expression with recurrence-free survival.

Breast cancers from three data sets were grouped based on if the expression of IGHM, IGHA2, IGJ, and IGHG1 mRNA was larger (red curves) or smaller (black curves) than the median expression. The figures display Kaplan–Meier curves for the indicated mRNAs (top titles) and cohorts (titles to the left) using recurrence-free survival as end point. The follow-up time on the x-axis is indicated in years. The p values were estimated with the log rank test. The y-axes indicate years after diagnosis and number at risk in the groups.

Fig. 6. Association of IGHA2 mRNA with prognosis upon adjustment for immune cell metagenes.

Multivariable Cox proportional hazards models of all cases and limited to basal-like case were performed for the SCAN-B cohort adjusting for different immune metagenes. The metagenes, indicated to the left in the figure, were generated using the microenvironment cell population Counter as for Fig. 1, as the mean log2 expression of mRNA-encoding cytokines specific for Th cells, or utilizing the log2 expression of FOXP3. The models were also adjusted for tumor size, lymph node status, patient age and chemotherapy treatment. For the model using all cases, adjustment was also done for PAM50 subgroup. The hazards ratios with 95% confidence interval for normalized log2 IGHA2 mRNA expression are shown in the figure.