An integrative microenvironment approach for follicular lymphoma: roles of inflammatory cell subsets and immune-response polymorphisms on disease clinical course

Clinicopathological features

A total of 237 FL patients were obtained (117 from UNICAMP and 120 from A. C. Camargo Cancer Center). The median age at diagnosis was of 57 years old (range: 19-94). As expected, there was a female predominance (133/237 or 56.1%), and most of the patients (189/237 or 79.7%) were diagnosed in advanced Ann Arbor stages (III or IV). The most common first-line treatment regimens employed were R-CHOP (107 patients) and R-CVP (65 patients). These and other characteristics are listed in Table 1.

Immunohistochemical assessment of cell components of the TME

Immunohistochemical (IHC) analysis of the TME components was possible in samples from 169/237 patients (71.3%), which were provided as a tissue microarray (TMA). The most frequent subpopulation was T-lymphocytes (median = 21.37% of pixels), with a predominance of cells expressing CD4 (median = 18.97%) over CD8 (median = 6.09%). The least frequent subpopulation was perforin-positive cells (median = 0.28% of pixels). Regarding the expression of cytokines, IL-17A had the strongest staining score (median H-score = 209.84), and IL-12A, the weakest (median H-score = 26.47). The antibody against TGFBR2 did not pass the quality control check in the tumor samples due to non-specific staining and thus, was not further analyzed. Figure 1 illustrates representative IHC stains before and after the pixel count.

The qualitative TME assessment included an evaluation of the meshworks of follicular dendritic cells (FDCs) and the infiltration patterns of TILs. Thirty-three patients (33/147 or 22.4%) were categorized as having no detectable FDCs (Group 1). Most of the tumors (64/147 or 43.6%) were classified as having poorly developed FDC meshworks (Group 2), and 34/147 patients (23.1%) were included in Group 3, which comprised predominantly well-developed meshworks. Finally, 16/147 patients (10.9%) had uniformly well-developed meshworks (Group 4) (representative IHC images are in Supplementary Figure 1).

An assessment of the patterns of the TILs revealed that the CD3+, CD4+, CD8+ and FOXP3+ cells presented a non-follicular infiltration pattern in most evaluable cases (89.3%, 87.2%, 90.2% and 69.6%, respectively). In contrast, for PD1+ and CD57+ cells, the follicular pattern was more prevalent (77.9% and 60.0% of evaluable cases) (Supplementary Table 1, Supplementary Figure 2). Moreover, significant differences in cell quantities were detected between patterns. There were higher quantities of CD3+ and CD8+ cells in cases that presented with a non-follicular pattern (p = 0.0001 and 0.03, respectively); whereas for PD1 and CD57, there were increased cell counts in patients with a follicular pattern (p = 0.02 and 0.01, respectively) (Figure 2).

The presence of a follicular pattern in the CD8+ cells was associated with bone marrow infiltration at diagnosis (p = 0.04). However, the patterns of the other TILs did not correlate with clinical features (Supplementary Table 2).

The quantifications (pixel counts) of the pan-markers (CD3 and CD68) revealed positive correlations with most other proteins. Conversely, the presence of NK cells (CD57+) correlated inversely with Ki-67, IL-12A, iNOS, and IL-10 expression. Furthermore, there was a negative correlation between the presence of CD8+ and granzyme B+ cells, indicating that not all CD8+ cells expressed cytotoxic markers. Finally, IL-17F correlated negatively with IL-12A, iNOS and granzyme B. The correlation matrix is shown in Figure 3.

The quantification of the TME cells had no impact on the clinical characteristics (data not shown). Regarding the expression of cytokines, higher IL-10 tumor levels were associated with a lower frequency of extranodal disease at diagnosis (p = 0.03). We also observed a trend of an association of high IL-10 levels with B-symptoms (p = 0.06) and with high-risk Follicular Lymphoma International Prognostic Index (FLIPI) (p = 0.05). The remaining cytokines were not associated with the clinical characteristics (Table 2).

Hierarchical clustering and principal component analysis of the TME components

To evaluate the TME proteins using a more realistic approach, we performed unsupervised hierarchical clustering of the immunohistochemical markers that were quantified using the Aperio system. Two clusters (1 and 2) resulted, with 100 and 32 patients, respectively (Figure 4A). Cluster 1 was enriched in cells that expressed CD8, perforin, CD57, PD1, and FOXP3 and exhibited higher Ki-67 counts than cluster 2. There was also higher expression of IL-17A, IL-17F and TGFBR1 for cluster 1. In contrast, cluster 2 had higher levels of granzyme B, IL-10, IL-12A, and iNOS (Figure 4B and 4C). No differences in macrophage infiltration (CD68, CD163) were seen (data not shown).

The clinical features of clusters 1 and 2 were similar with respect to the FLIPI index, B-symptoms, and bone marrow infiltration rates. The distributions of gender and age were similar as well. However, cluster 1 had a higher frequency of extranodal disease (p = 0.02) (Table 3).

Additionally, we performed an exploratory principal component analysis (PCA) using the IHC variables (both categorical and numerical data) to complement the strictly numerical data in the clusters. The final PCA model comprised 73.01% of the variance and was composed of 6 components (Table 4). The first component showed that elevated expression of TGFBR1, IL-12A, IL-2, and IL-10 was associated with a follicular pattern in FOXP3+ cells. The second component revealed a direct association between macrophages (CD68, CD163) and perforin expression. The third component aggregated the non-follicular proliferation patterns of CD3, CD4, and CD8 cells. The fourth component revealed an association of the follicular pattern of CD8 cells with elevated intratumoral IL-10, IL-17A, and IL-17F levels. For the fifth component, high CD68 and iNOS levels were associated with a high proliferative index (Ki-67). Finally, the sixth component demonstrated an opposing trend of high TGFβ and low IL-17A levels.

SNP genotyping and associations with clinicopathological features

A total of 159 FL patients (71 men and 88 women) were genotyped. There was a lack of Hardy-Weinberg Equilibrium (HWE) in 4 SNPs: IL12A rs568408, IL2 rs2069762, IL10 rs3024491, and TGFBR1 rs334348 (Supplementary Table 3).

Some SNPs were predictors of the clinical features at diagnosis. The genotypes GT+TT of IL2 rs6822844, compared with the GG genotype, were predictors for bulky disease at presentation (odds ratio [OR] = 2.55, 95% confidence interval [95% CI] = 1.10 – 5.91, p = 0.029). In addition, the presence of the CC genotype for IL10 rs1800872, in comparison with the CA+AA genotypes, was associated with bone marrow infiltration (OR = 2.07, 95% CI: 1.07–4.01, p = 0.030). Regarding IL10 rs1800890, the genotypes TA+AA were associated with the presence of B-symptoms when compared with the TT genotype (OR = 2.58, 95% CI = 1.29–5.15, p = 0.007). Finally, the AA genotype of TGFB1 rs6957, compared with the AG+GG genotypes, predicted the presence of extranodal disease (OR = 2.17, 95% CI = 1.02–4.58, p = 0.042).

Some of the SNPs were associated with the infiltration patterns of TILs in the FL biopsies (Table 5). For instance, the follicular pattern of the FOXP3+ cells could be predicted by the GA+AA genotypes of IL12A rs568408 (OR = 3.07, 95% CI = 1.10–8.56, p = 0.03). In addition, the SNP IL10 rs1800872 was associated with the follicular patterns of CD8 (CA+AA vs CC, OR = 8.87, 95% CI = 1.03–76.10, p = 0.04), PD1 (CA+AA vs CC, OR = 3.59, 95% CI = 1.17–11.0, p = 0.02) and CD57 (AA vs CC, OR = 8.30, 95% CI = 1.06–64.80, p = 0.04).

Other SNPs also associated with the quantification of TME proteins. First, a higher expression of iNOS was seen in the samples from FL patients carrying the “A” allele of IL12A rs755004 (p = 0.005, Figure 5A). Second, a higher percentage of granzyme B+ cells was observed in tumors from patients with the TT genotype of IL17F rs763780 (p = 0.04, Figure 5B). However, none of the SNPs altered the levels of the proteins that were encoded by each of their respective genes (data not shown).

Survival analyses

One-hundred seventy-two patients received R-CHOP or R-CVP as a first-line treatment and therefore were included in the survival analyses. The median follow-up times were of 6.05 years (all patients) and 6.32 years (living patients). At the end of the follow-up, 136 patients were alive (79.06%) and 36 (20.94%) had died.

We first evaluated the role of the main clinical features in event-free survival (EFS) and overall survival (OS) by estimating their hazard ratios (HRs) by univariate Cox regressions. As expected, the presence of high-risk FLIPI was associated with worse EFS (HR = 1.77, 95% CI: 1.08–2.88, p = 0.02) and OS (HR = 3.29, 95% CI: 1.64–6.59, p = 0.001) by univariate analysis. The same was observed for the presence of B-symptoms (HR = 2.39, 95% CI: 1.49–3.86, p < 0.001 for EFS and HR = 4.00, 95% CI: 2.02–7.93, p < 0.001 for OS). The presence of extranodal disease (excluding bone marrow) was also associated with worse EFS (HR = 1.86, 95% CI: 1.15–3.01, p < 0.01) and OS (HR = 2.48, 95% CI: 1.23–4.98, p < 0.01). Bone marrow infiltration was marginally associated with worse OS (HR = 1.89, 95% CI: 0.96–3.74, p = 0.06), but not with EFS (HR = 1.39, 95% CI: 0.87–2.22, p = 0.15). Finally, the administration of maintenance rituximab was associated with prolonged OS (HR = 0.40, 95% CI: 0.19–0.85, p = 0.01), but did not affect EFS (HR = 0.67, 95% CI: 0.41–1.09).

For the univariate analysis, the presence of follicular patterns for CD8 and CD57 was associated with worse EFS (p = 0.05 and 0.03, respectively); whereas, the follicular pattern of FOXP3 was associated with both worse EFS (p = 0.04) and OS (p = 0.01) (Kaplan-Meier curves on Figures 6A–6C and 7A). In addition, high immunohistochemical expression of IL-17F and the presence of FDC meshworks (Groups 2, 3 and 4) were associated with worse EFS (p = 0.01 for both markers, Kaplan-Meier curves on Figure 6D and 6E). Finally, a trend toward improved EFS (p = 0.06) and OS (p = 0.08) was seen in cases with high number of cells that expressed granzyme B (Figures 6F and 7B). After the multivariate analysis, the immunohistochemical variables that were independently associated with worse EFS included a follicular pattern of CD8+ cells (p = 0.001) and elevated expression of IL-17F (p = 0.04). Conversely, the presence of a follicular pattern of FOXP3+ cells was the only variable that was independently associated with worse OS (p = 0.02) (Table 6).

The prognostic roles of all SNPs and haplotypes were also evaluated. Linkage disequilibrium (LD) analyses allowed for the construction of haplotypes in IL10 and IL12A (Figure 8; frequencies of all the haplotypes in Supplementary Table 4). The AA-genotype of the polymorphism IL10 rs1800872 was slightly associated with worse OS by univariate analysis (p = 0.06, Figure 7C). In addition, the presence of a haplotype on IL12A (AA), involving the rs583911 and rs568408 “A” alleles, was associated with worse EFS and OS (p < 0.001 and p = 0.003, respectively; Kaplan-Meier curves in Figures 6G and 7D). After the multivariate analysis, the AA haplotype remained significant for EFS (p = 0.01) and was marginally associated with OS (p = 0.07) (Table 6).

Finally, we examined whether variables from the dimension-reducing analyses (PCA and clusters) altered the survival rates. In this setting, the presence of the third principal component was associated with prolonged EFS. In other words, patients with non-follicular infiltration patterns for CD3, CD4 and CD8 lymphocytes had better EFS (p = 0.02), and the statistical significance was maintained after the multivariate analysis (p = 0.01, Table 6). Complementarily, patients in cluster 1 had poorer EFS (p = 0.01) and a tendency toward worse OS (p = 0.06) than those in cluster 2 (Figures 6H and 7E). After adjusting for FLIPI, cluster 1 remained an independent predictor of worse EFS (p = 0.03) (Table 7).

Patients, controls and clinical data

The present retrospective study included an initial number of 237 cases, diagnosed with systemic FL (grades 1, 2, or 3A) between 1999 and 2016 at the Hematology and Hemotherapy Center – UNICAMP and A. C. Camargo Cancer Center.

The clinical data were collected from the medical files. FLIPI was used as the reference prognostic tool.

The study was approved by the Ethics Committees of both institutions (number 32177014.3.0000.5404), and all procedures were carried out according to the Declaration of Helsinki.

SNP selection and genotyping

The SNPs were studied using a “candidate gene” approach. SNPs located in immune response genes from the Th1, Th2, or Th17 axis, with previous evidence of active roles in lymphoma, cancer or immune response, were selected. A minor allele frequency of 5% was required. Sixteen SNPs were then selected, based on previous association studies (Table 8) [13, 18, 19, 29, 33–51].

DNA extraction and genotyping

DNA was extracted from the peripheral blood of patients and controls using the lithium chloride technique [52]. The DNA yield (ng/uL) and purity (260/280 and 260/230 ratios) assessments were performed using a spectrophotometer (NanoDrop^® 2000, ThermoFisher Scientific). The final concentration of all samples was set to 50 ng/uL.

SNP genotyping was performed on the Taqman^® OpenArray^® System (which is based on quantitative PCR reactions). Briefly, the DNA samples were pipetted together with Taqman^® Openarray^® Master Mix on 384-well plates. The mixture was transferred to genotyping plates by the Openarray^® Accufill^TM system. Thermocycling was performed over 40 cycles, and the visualization of the polymorphic alleles was achieved by VIC™ and FAM™ fluorophores. A single reaction allowed for the simultaneous detection of all 16 polymorphic variants.

Immunohistochemistry

The immunohistochemical reactions were performed using TMA slides to address the TME composition. A broad panel of TME markers was selected, which encompassed subpopulations of TILs, macrophages and FDCs (CD3, CD4, CD8, CD57, FOXP3, PD1, granzyme B, perforin, CD23, CD68, CD163, iNOS). In addition, for each gene in the SNP approach, the respective protein was detected by immunohistochemistry (IL-12A, IL-2, IL-10, TGFβ, TGFBR1, TGFBR2, IL-17A, IL-17F). Details on the clones and dilutions of each antibody are listed in Supplementary Table 5.

The TMA slides were subjected to antigen retrieval using citric acid/pH 6.0 or EDTA/pH 9.0 buffers. Endogenous peroxidase activity was blocked using a 3% hydrogen peroxide aqueous solution for 20 minutes. The antigen-antibody reaction was performed overnight at 4°C and was then amplified using a third-generation polymer that was tagged with anti-mouse and anti-rabbit immunoglobulins and horseradish peroxidase (HRP; Novolink Polymer Detection System, Leica Biosystems, Newcastle Upon Tyne, UK). The reactions were developed with the 3-3′-diaminobenzidine (DAB) chromogenic substrate (Sigma, D5637, St. Louis, MO, USA). Hence, the positive cell structures (membrane, cytoplasm or nuclei) presented as a brown color. Positive and negative controls were run for each batch. The negative control was obtained by omitting only the primary antibody used for the reaction.

IHC stain quantification was done automatically using the Aperio Scanscope XT device and by averaging both of the TMA core scores. The Positive Pixel Count algorithm was used to grade the pixels as negative (N), low positive (LP), positive (P) and high positive (HP) (Figure 1). The inputs for the algorithm were a hue value of 0.1, a hue width of 0.5 and a color saturation threshold of 0.1 (for most cores). In cases with nonspecific background, the color saturation threshold was increased to 0.15 to minimize false-positive detection. For the antibodies that stained specific TME populations (e.g. CD68, FOXP3, CD3), the fraction of all the positive pixels was considered as the score ($\text{Score}=\frac{LP+P+HP}{LP+P+HP+N}$). For the antibodies that heterogeneously stained both the tumor and microenvironment cells (all the cytokine antibodies), we calculated the H-score [53] using the following formula:

$$\text{HSCORE}=3\left(\frac{HP}{LP+P+HP+N}\right)+2\left(\frac{P}{LP+P+HP+N}\right)+1\left(\frac{LP}{LP+P+HP+N}\right)$$

For the TIL subpopulations, we also performed “qualitative analyses” by describing the predominant pattern of proliferation [5, 7, 8]. In this setting, there were 4 possible patterns as follows: “intrafollicular” (TILs homogeneously localized inside neoplastic follicles); “perifollicular” (TILs concentrated in the periphery of neoplastic follicles); “interfollicular” (TILs more visualized in the interfollicular spaces); and “diffuse” (TILs well distributed randomly inside and outside the follicles). For statistical purposes, “intrafollicular” and “perifollicular” were both considered “follicular” patterns, whereas “interfollicular” and “diffuse” were reclassified “non-follicular.”

For the CD23 stains, we used a “pattern” approach. Thus, the FDC categories were as follows: absent (Group 1); a minority of neoplastic follicles with disrupted meshworks (Group 2); the majority of neoplastic follicles with well- developed meshworks (Group 3); and uniformly well-developed meshworks (Group 4) [54].

Statistical analyses

We performed the chi-squared (χ²), Fisher’s exact test, Mann Whitney’s test, and Spearman’s correlation index to assess the bivariate associations. When necessary, numerical variables were dichotomized as “high” and “low,” based on the median levels.

In the SNP study, we tested, for each variant, the HWE, using the χ² goodness-of-fit test. Pairwise LD analyses were performed to estimate the haplotype formation using Haploview 4.2 (https://www.broadinstitute.org/haploview/haploview). LD was measured by the disequilibrium coefficient (D′), and LD significance was considered at a D′ ≥ 80%.

An exploratory PCA with Varimax rotation was used to address the interactions among the IHC scores. The PCA included both the numerical and categorical variables (pixel counting and patterns of tumor-infiltrating lymphocytes, respectively).

Agglomerative, unsupervised hierarchical clustering was also performed using the numerical IHC data (pixel counting) and using Euclidean distances and Ward’s clustering algorithm.

Survival analyses were performed only for the patients that received R-CHOP or R-CVP as a frontline therapy. The OS was defined as the time from diagnosis until death from any cause or the last follow-up. The EFS was defined as the time from diagnosis until death from disease, disease progression or the last follow-up. The survival curves were plotted using the Kaplan-Meier method and were compared using a log-rank test. We further performed Cox univariate regressions for the variables that influenced survival based on the Kaplan-Meier curves. Finally, a Cox multivariate model was proposed, including variables with a p value ≤ 0.10 in the univariate analysis. The genotyping and IHC data were evaluated separately in the survival models. The prognostic roles of both the PCA and clusters were also examined.

A final p value < 0.05 was considered statistically significant. When necessary, corrections for multiple comparisons were performed using the Benjamini-Hochberg method.

A CONSORT diagram (Supplementary Figure 3) summarizes the workflow of the study, as well as the number of patients for both the IHC and SNP analyses.