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

The study of the tumor microenvironment (TME) in follicular lymphoma (FL) has produced conflicting results due to assessment of limited TME subpopulations, and because of heterogeneous treatments among different cohorts. Also, important genetic determinants of immune response, such as single-nucleotide polymorphisms (SNPs), remain underexplored in this disease. We performed a detailed study of the TME in 169 FL biopsies using immunohistochemistry, encompassing lymphocytes, macrophages, and cytokines. We also genotyped 16 SNPs within key immune-response genes (IL12A, IL2, IL10, TGFB1, TGFBR1, TGFBR2, IL17A, and IL17F) in 159 patients. We tested associations between SNPs, clinicopathological features and TME composition, and proposed survival models in R-CHOP/R-CVP-treated patients. Presence of the IL12A rs568408 “A” allele associated with the follicular pattern of FOXP3+ cells. The IL12A AA haplotype included rs583911 and rs568408 and was an independent predictor of worse survival, together with the follicular patterns of T-cells (FOXP3+ and CD8+) and high IL-17F tumor levels. The patterns of CD3+, CD4+ and CD8+ cells, displayed as a principal component, also associated with survival. Hierarchical clustering of the TME proteins demonstrated a cluster that was associated with worse prognosis (tumors enriched in IL-17A, IL-17F, CD8, PD1, and Ki-67). The survival of FL patients who were treated in the rituximab era shows a strong dependence on TME signals, especially the T-cell infiltration patterns and IL-17F tumor levels. The presence of the AA haplotype of IL12A in the genome of FL patients is an additional prognostic factor that may modulate the composition of T-reg cells in this disease.


INTRODUCTION
Follicular lymphoma (FL) is the most common low-grade non-Hodgkin lymphoma (NHL) subtype and is characterized by an indolent clinical course and frequent relapses [1]. The progression of FL is influenced by stromal signals [2]. In this setting, the study of the tumor microenvironment (TME), which is composed of nonneoplastic cells and related molecules, has provided important insights on the signaling mechanisms and prognosis of FL. Several studies have described the role of tumor-infiltrating lymphocytes (TILs), macrophages, and inflammatory cytokines in the modulation of the clinical course of FL [3][4][5][6][7][8][9]. However, these results are often non-reproducible, due to differences in the cell quantification methods (e.g., visual estimation, manual counting, automated counting) and also due to the variation in treatment patterns [10,11].
Genes that encode inflammatory molecules, such as cytokines, also interfere in lymphoma biology, due to their intrinsic role in the development of lymphocytes [12]. In this setting, genomic variations, such as single-nucleotide polymorphisms (SNPs), especially variants that are associated with the modulation of immune functions, are relevant for the development and progression of lymphomas [12][13][14]. For instance, SNPs in key inflammatory genes, such as IL10 and IL2, are more consistently studied in large NHL cohorts and have been implicated in disease risk or in prognosis in the pre-rituximab era, including certain FL cohorts [12,13]. Cytokines that are encoded by other immune response genes, such as IL12A, TGFB1, IL17A, and IL17F, also regulate the balance of T-cell subsets and T-cell exhaustion in B-cell lymphomas. Nevertheless, the respective SNPs have been examined in few studies and in a non-integrated fashion with other components of the TME [14][15][16][17][18][19].
No study has yet evaluated the function of SNPs within immune genes in the TME composition of FL. This is biologically plausible, especially in the case of functional variants. In addition, there is a need for more specific prognostic factors in FL. The efficacy of rituximab (anti-CD20), which is associated with better treatment outcomes, may be influenced by components of the TME [9,10]. Therefore a more detailed investigation of the microenvironment in patients who receive anti-CD20 is warranted.
This study aimed to verify whether SNPs in immune response genes modify the TME composition and clinical features of FL. Another goal was to test the prognostic impact of these SNPs and the TME components in a cohort of patients who have been treated with rituximabcontaining regimens.

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: . 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.
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 panmarkers (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
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  The p values were obtained using Chi-squared tests ( * ) and Fisher's exact test ( ** ).   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) Oncotarget 3161 www.oncotarget.com 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).

DISCUSSION
We found that in FL patients uniformly treated with R-CHOP/R-CVP, the T-cell infiltration patterns, more importantly than quantification of these cells, independently associate with prognosis. In addition, we described a novel adverse and independent prognostic role of the AA haplotype in IL12A.
We reproduced previous successful evaluations of T-cell architectural patterns using TMAs [5,7,20] and found that the follicular pattern for FOXP3+ cells (T-reg cells) was independently associated with worse OS. Farinha et al. (2010) reported that this parameter shortened the time to transformation of FL in the pre-rituximab era [5]. We hereby confirm that this adverse impact also exists in patients receiving anti-CD20 therapy. In FL, FOXP3+ cells suppress the immune function of other T-cells, Oncotarget 3162 www.oncotarget.com impairing immunosurveillance [21]. Also, FL cells secrete CCL22, which induces the recruitment of T-reg cells [22].

Hence, it is conceivable that the concentration of T-reg cells near the FL cells (follicular pattern) is a histological representation of these immunosuppressive interactions.
Our findings also demonstrate that the follicular pattern of CD8+ cells was an independent factor for worse EFS, which is a novel finding for FL. In addition, one principal component encompassing non-follicular proliferations of T cell subsets (CD3+, CD4+ and CD8+), showed a favorable impact on EFS, which complements what was observed inside the neoplastic follicles. Indeed, T-cells in the interfollicular compartment seem to have preserved immune synapses and enhanced granzyme B production, possibly due to the involvement of the "border patrolling" effect [7,23,24]. In this setting, our results are plausible and corroborate the findings of Wahlin et al. (2010), associating non-follicular CD8+ cells with a good outcome in FL patients who have been treated in the prerituximab era [7].
Conversely, the presence of CD8+ lymphocytes in direct contact with FL cells (i.e. follicular proliferation pattern) promotes a defective immunological synapse due to the impairment of the CD8+ cell's F-actin cytoskeleton [25]. In addition, a subset of T-CD8+ cells in FL, located in the follicular compartment, may co-express molecules, such as ETV1, impairing immune function [26]. Therefore, it is plausible that CD8+ cells that are arranged in follicular patterns represent anergic subpopulations, which is reinforced by our finding associating the follicular patterns of CD8+ cells with worse outcomes.
Our initial survival analysis, as discussed above, showed that the T-cell patterns were associated with prognosis. However, none of the individual cell quantities significantly influenced survival. Notably, when we collectively evaluated these cells by hierarchical clustering, it became evident that the patients from cluster 1 had a poorer survival, independent of the FLIPI. Cluster 1 was enriched in cells with a high proliferation index (Ki-67) and in some T-cell subsets (CD8+, PD1+, FOXP3+). This cluster also presented with elevated expression of TGFBR1, IL-17A, and IL-17F. These results suggest that the quantities of the inflammatory cells, including T-cells, are important when they interact with each other and with cytokines, which reflects the complexity of TME in FL [27]. The mechanisms of these interactions should be explored further in future studies.
However, IL-17F was overexpressed in cluster 1 and was also found individually as an adverse prognostic factor for EFS. A previous report associated IL-17F expression with a bad clinical outcome of T-cell lymphomas, whereas for B-cell neoplasms, the scenario was unclear [14,15]. In this study, we described a novel adverse prognostic role for IL-17F in FL patients, mirroring what has been observed in T-cell lymphomas. The mechanism of this association has not been described, but it may reside in the downstream effects of IL-17F, such as angiogenic sprouting and prostaglandin production [14].     To complement the TME investigation in FL, we examined the genetics of the host. The AA haplotype of IL12A (constituted by rs583911 and rs568408 "A" alleles) was independently associated with worse EFS and was marginally associated with OS. Two previous studies did not find any association between rs568408 and outcome in FL. However, the cohorts in these studies were predominantly from the pre-rituximab era, limiting a direct comparison with our results [13,28]. The explanation for the prognostic role of the AA haplotype may reside in the association between the "A" allele of rs568408 and the follicular pattern of the distribution of the FOXP3+ cells. However, the exact mechanism of this association has not been elucidated. It is known that the rs568408 "A" allele may lower the activity of IL-12A by disrupting exonic splicing sites [29]. In this setting, a microenvironment with less active IL-12A may favor the outgrowth of fully immunosuppressive T-reg cells, which in FL, seems to be represented by such an arrangement as a follicular pattern, according to our results and others [5,30]. It is possible that other SNPs modulate the TME composition without altering the course of FL. For instance, the "A" allele of IL12A rs755004 was associated with increased intratumoral iNOS levels, which may be explained by the role of IL-12A on IFNγ secretion, with the consequent recruitment of fagocytes and oxidative burst [31]. We also found that IL17F rs763780 was associated with granzyme B levels, which may relate to the interference of Th17 cells on cytotoxic lymphocytes [32].
The most significant limitation of our study was its retrospective design. Therefore, the validation of our results using prospective cohorts is warranted. However, several strengths must be recognized. For example, our study had a relatively large sample size, and the patients were homogeneously treated with rituximab for our survival analyses, minimizing any therapy-related biases. We also performed computer-assisted image analyses for the quantification of the IHC markers, which has been proven superior and more reproducible in comparison to simple manual counting [11].
In this study, we demonstrated the adverse prognostic role of several TME elements in FL, considering both protein and genetic data. Importantly, the follicular pattern of the T-CD8+ cells, the high expression of IL-17F (individually and as a cluster member), and the AA haplotype of IL2A were all independently associated with worse prognosis. The latter variable was also associated with the pattern of FOXP3+ cell infiltration, which was validated in this study as a prognostic factor in the rituximab era. We believe that the simultaneous study of the tumor biopsies and immune response related SNPs in patients' peripheral blood allowed for a novel and more integrative approach that will provide new insights on follicular lymphoma biology.

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.

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™ 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.
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 (Score = 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 (χ 2 ), 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 χ 2 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 followup. 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.

Author contributions
Guilherme Assis Mendonça: designed the project, performed experiments, analyzed the data and wrote the manuscript. André Fattori: analyzed the data and wrote the manuscript. Rafael Rocha: analyzed the data and wrote the manuscript. Gustavo Lourenço: analyzed the data and wrote the manuscript. Márcia Torresan Delamain: analyzed the data. Suely Nonogaki: performed experiments and wrote the manuscript. Vladmir Lima: analyzed the data and wrote the manuscript. Gisele Colleoni: analyzed the data and wrote the manuscript. Cármino Souza: analyzed the data and wrote the manuscript. Fernando Soares: analyzed the data and wrote the manuscript. Carmen Lima: designed the project, analyzed the data and wrote the manuscript. José Vassallo: designed the project, analyzed the data and wrote the manuscript.