Association of body mass index and inflammatory dietary pattern with breast cancer pathologic and genomic immunophenotype in the nurses’ health study

The Nurses’ Health Study (NHS) was established in 1976 with the enrollment of 121,700 US female registered nurses aged 30–55 years, while the Nurses’ Health Study II (NHSII) consists of 116,429 US female registered nurses aged 25–42 years who were enrolled starting in 1989. Detailed descriptions of cohort procedures have been reported elsewhere [16]. Briefly, cohort members completed baseline questionnaires that provided medical histories and extensive information about demographic, lifestyle, reproductive, and dietary risk factors for breast cancer. Both cohorts have been followed biennially by mailed questionnaire to update information on exposure status and ascertain newly diagnosed diseases, including cancers. All women reporting incident diagnoses of breast cancer were asked for permission to review their medical records; 99% of reported cases were confirmed with pathology reports and medical record review. Tumor characteristics were extracted from pathology reports. Immunohistochemical evaluation of estrogen receptor (ER), progesterone receptor (PR), and human epidermal growth factor receptor 2 (HER2) expression was obtained from central review of breast tissue microarrays [16]. We analyzed 882 participants enrolled in the NHS and NHSII and diagnosed with invasive breast cancer, who had detailed exposure and follow-up data and gene expression data [9, 17]. The study protocol was approved by the institutional review boards of the Brigham and Women’s Hospital and Harvard T.H. Chan School of Public Health, and those of participating registries as required.

The empirical dietary inflammatory pattern (EDIP) score is a food-based dietary index for assessing dietary inflammatory potential, and it is based on circulating concentrations of C-reactive protein, interleukin-6, and tumor necrosis factor alpha receptor 2. EDIP is a weighted sum of 18 food groups from food frequency questionnaires (FFQs) derived in the NHS [14] and validated in NHSII and Health Professionals Follow-Up Study [13, 14]. The component food groups comprising the EDIP score that are positively associated with concentrations of inflammatory markers are: processed meat, red meat, organ meat, non-dark fish, vegetables other than green leafy vegetables/dark yellow vegetables, refined grains, high energy beverages (carbonated beverages with sugar, fruit drinks), low-energy beverages (low-energy cola/other carbonated beverages), and tomatoes. The component food groups that are inversely associated with concentrations of inflammatory markers are: beer, wine, tea, coffee, dark yellow vegetables (carrots, yellow squash, and sweet potatoes), green leafy vegetables, snacks, fruit juice, and pizza [14]. As previously described, the EDIP score in NHS/NSII ranges from − 3.34 to 2.81, with higher scores associated with higher concentrations of inflammatory markers [14]. Cumulative average EDIP uses all prior available FFQs (median 8, range 5–8 FFQs).

Weight and height were reported on study questionnaires, and self-reported weight and height were previously validated in the NHSII [18]. Height was obtained at enrollment, while weight (and other covariates) was updated every 2 years, starting in 1976 (NHS) or 1991 (NHSII). BMI (kg/m²) was calculated using weight from the participant’s last available questionnaire before breast cancer diagnosis (i.e., within 2–4 years of diagnosis). BMI change since age 18 was calculated as the last available pre-diagnosis BMI minus BMI at age 18. BMI change since age 18 provides a potential surrogate for adiposity relative to BMI alone, as significant rise in BMI from age 18 typically reflects increase in non-lean body mass [19]. Prior studies suggest adult weight gain captures dynamic pattern of weight trajectory throughout adult life reflecting a time-integrated metric, adults gain weight mostly through accumulating fat mass and detrimental fat distribution, and adult weight gain is a simpler and more intuitive concept to the general public than BMI [19]. As a sensitivity analysis, the last available pre-diagnosis BMI was converted to a categorical variable with three levels (underweight/normal weight, overweight, and obese) as determined by the standard BMI thresholds for women (< 25, 25–29.9, > 30). Very few participants had a BMI category of “underweight” (N = 11/882, 1.2%) so the “underweight” and “normal weight” categories were combined for this analysis. Further sensitivity analyses assessed pre-diagnosis BMI as a continuous predictor and BMI at age 18 as a continuous predictor. For physical activity, we used cumulative average metabolic equivalent task (MET) hours per week to capture long-term activity. Other covariates such as race, age at first birth, parity, age at diagnosis, year of diagnosis, menopausal status, recent postmenopausal hormone therapy (HT), and smoking [20] were retrieved from baseline, subsequent, or most recent NHS/NHSII questionnaires. When available, the most recent pre-diagnosis questionnaire was prioritized.

Tumor microarrays (TMA) were constructed from archived formalin-fixed paraffin-embedded (FFPE) breast tumor blocks at the Specialized Histopathology Core of the Dana Farber/Harvard Cancer Center with three 0.6 mm core biopsies taken from the tissue blocks for each tumor [16]. IHC was conducted on 5 μm paraffin sections of the TMA blocks for CD8 (Dako 7103, clone C8/144B), CD4 (Dako 7310, clone 4B12), CD20 (Dako 0755, clone L26), and CD163 (Vector Labs VP-C374, clone 10D6) [21]. Staining positivity was determined using an automated computational image analysis system (Definiens Tissue Studio software, Munich, Germany). The percent of positive cells for each antibody was determined as the number of cells positive over the total number of cells (nuclei) for each of the three cores from each patient. The final mean percent positivity was the average of the three cores.

For NHS, RNA was extracted from three cores of 1 mm or 1.5 mm diameter taken from FFPE blocks of tumor tissue as described previously [22]. RNA expression was determined using Affymetrix Glue Grant Human Transcriptome Array [23]. Gene expression data were normalized and summarized using robust multiarray average (Affymetrix Power Tools v1.18.0), correcting for batch effects using a surrogate variable analysis R package, ComBat, and excluding low expressing genes (< 25th percentile) [9, 17]. The final dataset consisted of 15,369 annotated probesets of coding and non-coding RNAs. All microarray and annotation data are available at the National Center for Biotechnology Information Gene Expression Omnibus (accession number: GSE115577). For The Cancer Genome Atlas (TCGA) [24], gene expression data were obtained from XENAbrowser (version 2015-02-24) for breast cancer samples with paired tumor and normal expression data.

Published immune signature scores were derived as the median value of the genes in each gene set from gene expression data for each patient’s tumor, a validated approach for analyses of disparately derived signatures [25]. This included 9 gene expression signatures demonstrated to correlate with infiltrating immune cells [26‐30], including the GeparSixto (GSAct) 6-gene signature (CXCL9, CCL5, CD8A, CD80, CXCL13, CR2), which was derived based on association with tumor infiltrating lymphocytes in the neoadjuvant chemotherapy breast cancer clinical trial GeparSixto [26]. In addition, we calculated the immune score for each of the 96 gene sets in the Immune Response In Silico (IRIS) signatures [31], a compendium of immune cell-specific gene sets derived to reflect immune cell activation/exhaustion states, available in ImmuneSigDB [32].

For each of the four immune IHC markers of interest (CD4, CD8, CD20, and CD163), we calculated total mean percent positivity as described above. As none of the published immune signature scores individually demonstrated robust association with IHC, we evaluated linear models (four models, one for each immune cell subset) composed of all 105 published immune signatures available as predictors and the mean percent IHC positivity for each subset as the outcome. We then performed lasso reduction on these models, resulting in linear models of signatures that are significantly associated with the mean percent of IHC positive immune cells for each subset. We used the models to calculate an immune cell-specific score for each patient’s tumor, derived as a weighted arithmetic sum of the linear model for each patient’s gene expression data. The predicted scores represent the predicted amount of that particular IHC marker, with higher score corresponding to higher IHC. These scores were the main outcomes of interest for the present analyses (analyses schema in Fig. 1).

We evaluated four immune scores based on the lasso reduction of IRIS gene expression signatures that were associated with IHC for four immune cell-specific markers: CD4+, CD8+, CD20+, and CD163+. We assessed the correlation of each of the immune scores with the other three scores. As an established comparison, we also assessed the GSAct. Correlations were measured using Spearman rho correlation coefficients and corresponding p values. Bivariable and multivariable linear models were used to compare the CD4+, CD8+, CD20+, and CD163+ immune scores (outcomes) with BMI change since age 18, physical activity, and EDIP score (predictors). Covariates included in the multivariable analyses are detailed above. All analyses and figure generation were performed using R version 3.5.1.

A total of 882 subjects enrolled in NHS or NHSII with available gene expression data derived from tumor samples were included in the study (Table 1; Fig. 2A). Overall, most patients in the cohort were white (95.5%; 842/882), stage I–II at diagnosis (90.6%; 799/882), and the majority of patients had ER+/HER2−negative breast cancer, followed by ER±/HER2+ and TNBC (60.2%, 27.4%, and 12.4%, respectively) (Additional file 1: Table S1). Of the gene expression cohort, 623 also had gene expression data derived from adjacent normal breast tissue and 262 subjects had IHC data for at least one immune marker (Figs. 1, 2A). Participants who had tumor gene expression data plus IHC data (the “IHC training cohort”), when compared to those with expression data alone (the “expression only evaluation cohort”), were exclusively from NHSII. Thus, based on cohort differences, the IHC training cohort participants were more likely to be diagnosed after 2000 (58.0% vs. 45.2%) and less than age 60 at diagnosis (100% vs. 31.0%) (Table 1). Compared to the expression only cohort, the IHC cohort had a greater proportion of ER+/HER2− breast cancers (68.3% vs. 50.5%) and more patients with normal BMI (57.6% vs. 46.6%). There was no significant difference in stage at diagnosis, cumulative average physical activity, or cumulative average EDIP score by cohort (Table 1).

To understand the landscape of immune gene expression, in the full cohort (n = 882 participants) we first evaluated 105 previously published, well-established immune gene expression signatures representing general immune infiltration as well as specific immune subsets (Additional file 2: Table S2). Among these, we visualized nine general immune signatures shown to be associated with TILs [26‐30] (Fig. 2B). These signatures indicated that most tumors had relatively low expression of immune-related genes, while a small subset (~ 10–15%) showed high immune gene expression.

Among tumors with paired IHC and gene expression (IHC training cohort), we evaluated the association of each of the 105 published immune signatures with IHC staining for CD8 (n = 250), CD4 (n = 258), CD20 (n = 252), and CD163 (n = 253) (Additional file 3: Table S3). Correlations were modest, with most Spearman correlation coefficients < 0.30 and none greater than 0.40. Among all signatures, the GeparSixto “immune activation signature” (GSAct), which was derived based on association with TILs in the GeparSixto clinical trial [26] and was recently confirmed to have a robust association with response to neoadjuvant chemotherapy in early-stage breast cancer [33], was among the strongest associations across all immune IHC markers (Additional file 3: Table S3). This suggested that, as designed, GSAct was a representative signature of TILs in breast tumors; thus, moving forward in analyses, we focused on GSAct. As anticipated, GSAct was highest in TNBC relative to other breast cancer receptor subtypes as well as higher in basal-like and HER2-enriched molecular intrinsic subtypes (Additional file 4: Figure S1A, B).

As an initial validation, we evaluated the association of GSAct with CD8+ IHC among 250 patients with paired CD8+ T cell immunohistochemistry and gene expression data. Overall, CD8+ quantification corresponded to visualized CD8+ cell amount (example images in Fig. 2C); this supports recent work from our group showing overall good agreement between computational evaluation of IHC with expert pathology review (Roberts, Tamimi, et al.; manuscript under review).

When stratified into quartiles by reported total CD8 IHC quantification, there was a strong association with GSAct (ANOVA p = 1.9e−7; Fig. 2D; CD8 IHC continuous; Additional file 4: Figure S1C-D). The overall correlation with CD8 IHC and GSAct as continuous variables remained for epithelial, stromal, and total CD8 positivity (range Pearson’s r = 0.45–0.47; Additional file 4: Figure S1E).

The NHS cohort offers a unique opportunity to interrogate gene expression data from tumor and adjacent normal tissue. We focused our comparisons on GSAct and CD8A as a representative single gene for general CD8+ T cell infiltration. Among the 623 patients with paired tumor–normal samples, we observed a positive correlation between tumor and adjacent normal tissue for GSAct (Pearson’s r = 0.31, p < 0.001) and CD8A single gene (Pearson’s r = 0.46, p < 0.0001) (Fig. 3A). These results were validated in the independent TCGA breast cancer cohort, with 113 pairs of tumor and adjacent normal breast transcriptome data (Pearson’s r = 0.45, p < 0.0001; Pearson’s r = 0.34, p = 0.0002 for GSAct and CD8A, respectively; Fig. 3B). Of the adjacent normal tissue samples in NHS, 297 were classified as a PAM50 subtype other than “normal-like” [17]. As a sensitivity analysis, we removed these 297 samples, which could be a result of contamination of normal samples, and repeated all analyses (Additional file 4: Figure S2A-B), with no meaningful change in results. As an exploratory analysis, we assessed differences in paired tumor–normal samples among the NHS and TCGA cohorts using a paired Wilcoxon test. Overall, in the NHS cohort we found that GSAct and CD8A each differed significantly between tumor and normal samples (p < 0.0001 and 0.002, respectively). However, there did not appear to be a consistent trend when comparing tumor/normal by receptor subtype or PAM50 subtype for either GSAct or CD8A (Additional file 4: Figure S2C-F).

We initially attempted to derive novel immune signatures based on individual genes’ contributions; however, we found no significant discrimination relative to published immune expression signatures (Additional file 3: Table S3). Instead, we pursued an approach using all 105 published immune gene expression signatures and lasso reduction using the IHC training cohort. Eighteen signatures remained in the CD4+ lasso regression model, 24 remained in the CD8+ model, 24 remained in the CD20+ model, and 12 remained in the CD163+ model (Additional file 1: Tables S4–7).

We used the final models to calculate a patient-specific score for each immune subset, derived as a weighted arithmetic sum of the linear model for each patient’s gene expression data, and plotted the score against that patient’s mean IHC percent positivity for each marker to assess the fit of our model (Fig. 4A–D). The Spearman’s rho correlation statistic ranged from 0.42 to 0.54 with p values all < 0.0001. We compared the performance of our four immune cell-specific lasso model scores with the 105 individual immune expression scores (9 immune signatures plus 96 IRIS gene sets; Additional file 3: Table S3). As anticipated, for each IHC marker our novel immune scores demonstrated the strongest correlation with IHC quantification. To gain a descriptive understanding, we visualized all 882 tumors using our four novel IHC-based immune expression scores along with receptor and PAM50 subtype (Fig. 4E) and assessed the correlation between each immune score with each other (Additional file 4: Figure S3A). All immune scores were significantly, positively correlated with each other (Spearman’s rho ranging from 0.57 to 0.81). CD8+ with CD20+ scores had the highest correlation (Pearson’s r = 0.81, p < 0.0001). We further assessed these correlations by IHC subtype and found similar correlations (Additional file 4: Figure S3B-D). Similarly, CD8+ and CD20+ scores were most highly correlated with each other among all IHC subtypes (Spearman’s rho: 0.79 to 0.84, all p < 0.0001).

Using these novel IHC-validated immune expression scores, we evaluated the association with the immune microenvironment and patient features – BMI change since age 18, physical activity, and EDIP score in the expression-only evaluation cohort (n = 620 participants), omitting the IHC-only cohort as the dataset on which the signatures were derived. In bivariable and multivariable analyses, there was no significant association between cumulative average physical activity (MET hours/week) and GSAct, CD8+ score, CD4+ score, CD20+ score, or CD163+ score (Tables 2, 3). Among the IHC training cohort, there was no association between BMI change since age 18 and any immune score. Among the expression only evaluation cohort, CD4+ score, and CD163+ score were positively associated with BMI change since age 18 in bivariable (p = 0.003, and p = 0.04, respectively) and multivariable (p = 0.009, and p = 0.04, respectively) analyses. In other words, for each 10 unit (kg/m²) increase in BMI, the percentage of cells positive for CD4 and CD163 increased 1.6% and 1.4%, respectively. In the multivariable analysis only, BMI change since age 18 was positively associated with CD8+ score (p = 0.05). Cumulative average EDIP score was not significantly associated with CD8+, CD4+, CD20+, or CD163+ scores in either bivariable or multivariable analyses among either cohort.

We performed sensitivity analyses of pre-diagnosis BMI (one two-year follow-up cycle prior to diagnosis) as a categorical predictor, pre-diagnosis BMI as a continuous predictor, and BMI at age 18 as a continuous predictor (Additional file 1: Tables S8–10). In multivariable models, categorical BMI was not associated with CD8+ or CD163+ score in either cohort, but obesity, compared to normal BMI, was positively associated with CD4+ score (p = 0.02) among the expression only cohort. Among the IHC cohort, continuous pre-diagnosis BMI and BMI at age 18 were not associated with any immune score. Among the expression only cohort, continuous pre-diagnosis BMI was positively associated with GSAct (p = 0.05), CD8+ (p = 0.02), CD4+ (p = 0.003), and CD163+ (p = 0.03); however, BMI at age 18 was not associated with any immune score.