Introduction
Follicular lymphoma (FL) is the second most frequent non-Hodgkin lymphoma and is characterized by a generally indolent clinical evolution. While asymptomatic and low-tumor burden patients can initially be managed by observation, symptomatic patients and patients with high-tumor burden disease are typically managed at diagnosis with immunochemotherapy (IC) with or without maintenance therapy with an anti-CD20 monoclonal antibody as standard of care [1]. IC-treated patients who do not have disease relapse, progression, transformation, or re-treatment within the first 24 months after diagnosis can expect to have a subsequent mortality risk equivalent to the general population [2].
In contrast, patients who do not achieve event-free status (EFS) at 24 months (EFS24) have a very poor outcome [3].
There is an intense interest to identify FL patients at diagnosis who are at risk for poor outcomes with the development of prognostic scores that incorporate tumor-based biological variables such as the PRIMA 23-gene expression profiling (GEP) signature [4], tumor mutations (m7-FLIPI) [5], or other genomic alterations [6]. The prognostic role of host germline genetics has been less studied [7] and has mainly focused on candidate genes or pathways such as FcG receptors (FCγR) [8], complement [9], cytokines [10], and DNA repair and one-carbon metabolism [11]. In the only genome-wide association study (GWAS) of FL outcomes in 586 patients from 3 observational studies from the pre-rituximab era [12], a suggestive locus at 17q24 near ABC transporter ABCA10 and ABCA6 genes (P = 5.24 × 10−8) was identified.
Unlike GWAS for cancer etiology, where a large number of case-control studies are often available, GWAS for cancer prognosis have been hampered by fewer studies, smaller sample sizes of individual studies (particularly the number of events), heterogeneity of available studies, and a lack of data on clinical characteristics, treatment(s) and key clinical outcomes (e.g., EFS). To address these issues, we conducted a GWAS to identify loci associated with prognosis in FL patients treated with first-line IC (N = 1054) using three randomized phase 3 clinical trials and a prospective observational cohort of IC-treated patients. We conducted both a standard GWAS that pooled the four cohorts, as well as a complementary statistical approach inspired from cross-validation designs.
SNPs prioritized based on statistical significance and their ability to replicate across our “leave-one cohort out” (LOCO) design were further evaluated in a polygenic risk score (PRS) model and for in silico functional annotation (Fig. 1).Fig. 1Overview of the follicular lymphoma GWAS.
Image was created with BioRender.com.
Methods
Study design and population
We conducted a GWAS to identify germline genetic predictors of EFS as well as EFS24 after treatment initiation in FL patients treated with first-line IC. Patients (N = 1054) were from prospective trials of the Lymphoma Study Association, the Fondazione Italiana Linfomi, and the Molecular Epidemiology Resource (MER), a prospective observational cohort from the University of Iowa and Mayo Clinic. The clinical trial cohorts included patients with available GWAS data from the PRIMA [13] (N = 391), RELEVANCE [14] control arm (N = 205), and FOLL05 [15] (N = 199) trials.
The MER cohort [16] consisted of 259 FL patients prospectively enrolled in the MER. All diagnoses were confirmed by specialized hematopathologists and excluded grade 3b. Further details of the cohorts and treatments are provided in Supplementary Information.
Genotyping and imputation
Genotyping was performed using the Illumina Human Core exome (Illumina, San Diego, CA, USA) at the Centre National de Genetique (Evry, France) for PRIMA, RELEVANCE, and FOLL05. The MER was genotyped at the NCI Cancer Genomic Research Laboratory as part of two InterLymph initiatives, the first using the Illumina OmniExpress (N = 154) and the second using the Illumina Global Screening Array (N = 105) [17]. All genotyping used standard quality control measures.
Population substructure across the cohorts was evaluated using EIGENSTRAT [18], and plots of the first four principal components analysis (PCA) are provided in Supplementary Fig. 1; no outliers required removal. To harmonize genotypes that were obtained from different arrays and to include additional genotypes across the pooled GWAS analysis, we performed an imputation with Minimac4, using 1000G as the reference panel.
The quality control and filtering steps are summarized in Supplementary Table 1. The final GWAS analysis was performed on 1,287,537 high-quality SNPs common across the cohorts with a minor allele frequency (MAF) > 0.05.
Statistical analyses
GWAS of EFS and EFS24
EFS was the primary endpoint, and was defined as the time from the date of study entry (FOLL05, PRIMA, RELEVANCE) or the date of diagnosis for the MER to the date of disease progression/relapse, transformation, re-treatment, or death from any cause; if none of these occurred, we used the date of last update in the study [16, 19–21]. The secondary endpoints were the evaluation of SNPs associated with these events within or after the first 24 months after treatment initiation (or diagnosis), corresponding to EFS24 and OS. Association testing between SNPs with EFS and OS was conducted assuming a log-additive genetic model.
We calculated the hazard ratios (HR) and 95% confidence intervals (CIs) using Cox regression analysis. All models were adjusted for sex, cohort, four eigenvectors from the PCA, Follicular Lymphoma International Prognostic Index (FLIPI), treatment (R-CHOP [rituximab, doxorubicin, cyclophosphamide, vincristine, prednisolone]/R-FM [rituximab, fludarabine, mitoxantrone] vs. R-CVP [rituximab, cyclophosphamide, vincristine, prednisolone]), and use of rituximab maintenance (RM).
For FLIPI adjustment, we used two categories (0–2 and 3–5) because we did not have all the individual FLIPI variables in FOLL05. As RM began after induction IC, this variable was modeled as a time-dependent covariate with a cut-off of 6 months after the initial date of study entry (FOLL05, PRIMA, RELEVANCE) or the date of diagnosis for the MER. For treatment adjustment, we grouped R-CHOP and R-FM based on the results of the FOLL05 trial showing similar outcomes between the two arms but a better outcome than the R-CVP arm [15].
Finally, to identify SNPs associated with events within the first 24 months (EFS24) versus after 24 months, we used a piecewise Cox model [22].
Statistical model for GWAS
We applied two approaches for statistical analysis to better account for the heterogeneity of the four cohorts and address the robustness of the results (i.e., replication across cohorts) of genetic variants associated with the outcome. We first conducted a standard GWAS that pooled the four cohorts to obtain the distribution of P values of SNPs (observed and imputed) associated with EFS; P < 5 × 10−8 was defined as genome-wide significance (GWS) and P < 1 × 10−5 as suggestive. The second step consisted of a method inspired from cross-validation statistical approaches.
The data were iteratively analyzed in training and testing blocks by LOCO in each round, such that a LOCO-test set consisted of three cohorts, with the cohort left out used as the LOCO-test set (Fig. 2 and Supplementary Fig. 2).
The process was repeated such that each cohort was used iteratively, defining LOCO rounds #1 to #4. For each LOCO round, each SNP was scored +1 if and only if PTrain-set < 10−4 and PTest-set < 0.1 and there was a consistent HR direction in both the train and test sets (i.e., no HR flip); the SNP was scored 0 in all other cases. These P value thresholds for LOCO were set to balance both type 1 and type 2 errors, considering the available cohort sizes and based on asymptotic power calculations (Supplementary Table 2).
For each SNP, the total LOCO score was obtained by summing the scores obtained in the four LOCO rounds, ranging from 0/4 to 4/4. LocusZoom [http://locuszoom.org] was used to plot results for each sentinel SNP and associated region identified from the GWAS of EFS and EFS24.Fig. 2LOCO model description and application in this study.The dataset was split into a train set with three cohorts and a test set corresponding to one cohort, defining one LOCO round.
The process was repeated such that each cohort was used iteratively for the train and test phases, for a total of four LOCO rounds (#1 to #4). In each LOCO round, a SNP was scored +1 if and only if the P-train-set <10−4 and P-test-set <0.1 and there was a consistent HR direction in both the train and test sets (i.e., no HR flip); else the SNP was scored 0. For each SNP, the total LOCO score was obtained by summing the scores obtained in the four LOCO rounds, ranging from 0/4 to 4/4.
Of note, the P value thresholds for LOCO were set to balance both type 1 and type 2 errors, considering the available cohort sizes and based on asymptotic power calculations (see Supplementary Table 2).
Multi-SNP risk score for combined EFS and EFS24
Using the top independent SNPs from the GWAS, we generated a PRS corresponding to a weighted average of the risk alleles, the weights being the log of the adjusted HR for EFS and EFS24. Kaplan–Meier curves were generated to depict the multi-SNP score based on PRS quantiles, overall and by FLIPI subgroups.
Overall survival (OS) and histological transformation (HT)
In exploratory analyses, we tested our SNPs of interest from the EFS/EFS24 analyses for their association with OS. To investigate the association of HT with the top SNPs, HT progression without transformation and death (all defined at first event) were modeled as competing risks [23]. These analyses were performed only for PRIMA and MER cohorts in which HT was specifically documented as part of previous studies [19, 24, 25].
Multi-SNP risk score in cohorts with other first-line treatments
As exploratory analyses, the PRS was also tested for EFS for patients treated with the combination of rituximab and lenalidomide in the experimental arm of the RELEVANCE trial [14, 20]; patients from MER treated by rituximab in combination with bendamustine; and patients in the MER treated with rituximab monotherapy [16].
Multivariable survival model
Multi-SNP score based on PRS of the top SNPs was then integrated in a multivariable Cox model with FLIPI, sex, treatment, RM, and cohort. As all individual variables for FLIPI were available for PRIMA, MER, and RELEVANCE (0–5 categorization) but not in FOLL05 (0–2 and 3–5 categorization), the multivariable model used a “saturated” FLIPI to integrate all available FLIPI information. Finally, we used the Harrell index corrected for optimism using cross-validation to compare the prognostic ability of survival models [26].
Post-GWAS analyses
Correlation with RNA sequencing (RNAseq)
RNAseq data from FFPE diagnostic tumors of 222 patients from the RELEVANCE cohort, along with paired GWAS data, were available. Detailed methods for RNAseq have been published [27]. We used the GWAS and RNAseq data to assess the correlation between the genotypes of top SNPs and expression of transcripts (eQTL) of coding genes within 500 kb from each sentinel SNP using the positions from GenCode v24 (Supplementary Information).
Correlation with assay for transposase-accessible chromatin with high-throughput sequencing (ATACseq)
ATACseq data from 56 frozen diagnostic tumors of patients included in the RELEVANCE cohort were available for correlation between the different genotypes of top SNPs and open chromatin region (OCR) profiles. Frozen tumor tissue was processed according to the previously published FAST-ATAC and OMNI-ATACseq protocols [28] (Supplementary Information).
Functional annotation
We used in silico approaches to better understand the potential functional impacts of top SNPs associated with EFS and EFS24. Evidence of functionality was assessed using HaploReg [29] [https://pubs.broadinstitute.org/mammals/haploreg/haploreg.php], RegulomeDB [30] [https://regulomedb.org/regulome-search/] with specific research on “B-cells tissue” regarding the origin of FL and FORGEdb score [31] [https://forge2.altiusinstitute.org/files/forgedb.html].
Results
Study patients
Overall, the median age at FL diagnosis was 57 years (range, 22–87 years), 55.7% were male, and 40.9% of patients had FLIPI 3–5 score. The clinical characteristics of patients, type of IC, and use of RM by study cohort are summarized in Table 1. With a median follow-up of 99.5 months (95% CI: 95.4–107.4) in the pooled cohort, the 6-year EFS rate was 56.5% (95% CI: 53.6–59.6) and the 6-year OS rate was 89.5% (95% CI: 87.7–91.4), with results similar across individual cohorts (Table 1 and Supplementary Fig.
3). EFS was superior for FLIPI 0–2 vs 3–5 in all cohorts and with RM in PRIMA and MER (Supplementary Fig. 4).
Patients who failed to achieve EFS24 had a much worse subsequent OS compared to patients achieving EFS24 (Supplementary Fig. 5).Table 1Demographics, clinical characteristics, and outcomes of patients in the FL GWAS.CharacteristicsFOLL05 n = 199MER n = 259PRIMA n = 391RELEVANCE n = 205All n = 1054Recruitment period2006–20102002–20152004–20072011–2014AgeMean (years)55.357.256.459.7157.0Range (years)31–7523–8622–8733–8322–87Sex,n(%)Female84 (42.2)109 (42.1)167 (42.7)107 (52.2)467 (44.3)Male115 (57.8)150 (57.9)224 (57.3)98 (47.8)587 (55.7)FL grade,n(%)1 or 2161 (80.9)160 (61.7)333 (85.2)176 (85.9)830 (78.8)3a22 (11.1)99 (38.2)22 (5.6)11 (5.4)154 (14.6)Unspecified 1-3a16 (8)0 (0)36 (9.2)18 (8.8)70 (6.6)Ann Arbor stage,n(%)I0 (0.0)24 (9.3)5 (1.3)0 (0)29 (2.8)II16 (8.0)41 (16.0)27 (6.9)16 (7.8)100 (9.5)III60 (30.2)73 (28.4)77 (19.7)38 (18.5)248 (23.6)IV123 (61.8)119 (46.3)282 (72.1)151 (73.7)675 (64.2)Missing-2--2FLIPI score,n(%)a0–2119 (59.8)177 (68.3)229 (58.7)97 (47.3)622 (59.1)3–580 (40.2)82 (31.7)161 (41.3)108 (52.7)431 (40.9)Missing--1-1Induction therapy,n(%)R-CHOP63 (31.7)163 (62.9)308 (78.8)192 (93.7)726 (68.9)R-CVP65 (32.7)96 (37.1)81 (20.7)13 (6.3)255 (24.2)R-FM71 (35.7)0 (0)2 (0.5)0 (0)73 (6.9)R maintenance,n(%)Yes6 (3.0)53 (20.5)173 (44.2)205 (100)437 (41.5)No193 (97.0)206 (79.5)218 (55.8)0 (0)617 (58.5)OutcomesMedian FU, months83.7143.5119.281.099.595% CI81.5–86.7142.9–155.1117.7–121.578.2–83.195.4–107.4Number of events for EFS9513521089529Number of deaths for OS246672251872-year EFS rate73.876.177.280.877.095% CI67.9–80.271.0–81.473.1–81.575.5–86.474.5–79.66-year EFS rate54.259.653.859.456.595% CI47.6–61.753.9–66.049.0–59.152.9–66.853.6–59.66-year OS rate91.888.089.689.189.595% CI88–95.584.1–92.086.6–92.884.9–93.587.7–91.4FL follicular lymphoma, FLIPI Follicular Lymphoma International Prognostic Index, R-CHOP rituximab, doxorubicin, cyclophosphamide, vincristine, prednisone, R-CVP rituximab, cyclophosphamide, vincristine, prednisone, R-FM rituximab, fludarabine, mitoxantrone, R rituximab, FU follow-up, CI confidence interval, EFS event-free survival, OS overall survival.aFLIPI is measured at the date of diagnosis in MER and at enrollment for the three prospective trials (FOLL05, PRIMA, RELEVANCE).
SNP discovery and internal validation
After applying rigorous quality control filters and imputation, data on 1,287,537 SNPs with a MAF > 5% in the pooled cohort of 1,054 patients were available for analysis (Supplementary Table 1). We conducted both a standard GWAS analysis of the pooled cohorts and then overlaid those results with the results from an internal cross-validation design (LOCO, Fig. 2 and Supplementary Fig.
2) to identify the most reproducible genetic variants associated with EFS and EFS24.
SNPs associated with EFS
Examination of the quantile-quantile plot (Supplementary Fig. 6) showed minimal evidence for population substructure, overall (λ = 1.01). The Manhattan plot overlayed with the LOCO results (Fig.
3a) revealed suggestive associations but no SNPs meeting the GWS threshold, with strongest association signal for EFS at 3q27.3 (rs72625024, P = 9.37 × 10−8) with a LOCO score of 4/4 reflecting the robustness of the association with EFS across the four cohorts (HRs ranging from 0.54 in MER to 0.71 in RELEVANCE) and from the pooled analysis with an HR = 0.62 (Table 2 and Supplementary Table 3). This SNP is localized in HRG/FETUB antisense RNA 1 affiliated to a long-non-coding RNA (lncRNA) near the histidine rich glycoprotein (HRG) gene and FETUB (fetuin B) gene (Fig. 4a).
This was followed by loci at 3q21.1 (rs1440146, P = 2.24 × 10−6) with a LOCO score of 1/4 and an HR = 1.35 from the pooled analysis; 9q21.31 (rs10867429, P = 5.67 × 10−6) with a LOCO score of 0/4; and 8p23.2 (rs6558733, P = 7.44 × 10−6) with a LOCO score of 2/4 and an HR = 1.50 from the pooled analysis. The 3q21.1 SNP is in an intronic region of the Kalirin RhoGEF kinase (KALRN) gene (Fig. 4b), while the 8p23 SNP is in an intronic region of the CUB and Sushi multiple domains 1 (CSMD1) gene (Fig.
4c).Fig. 3Manhattan plot for genome-wide association study from the pooled analysis overlayed with the LOCO results.Data shown are from the combined meta-analysis, including 1054 patients. The chromosomes are ordered on the x-axis; the y-axis shows the −log10(P) values for the association tests with event-free survival for the whole period of follow-up (a), the period within the first 24 months (b), and the period after the first 24 months (c).Fig.
4Analysis of follicular lymphoma GWAS for event-free survival identifies variants near HRG/FETUB, KALRN, CSMD1, TCL6, WSCD2/CMKLR1, and PPIAP1 genes.Regional association and linkage disequilibrium (LD) plots at HRG/FETUB (a rs72625024), KALRN (b rs1440146), CSMD1 (c rs6558733) for event-free survival for whole period of follow-up and at TCL6 (d rs114695031), WSCD2/CMKLR1 (e rs12824955), KALRN (f rs333302) and PPIAP1 (g rs118086521) for event-free survival for the first 24 months of follow-up (EFS24). The y-axis represents the −log10(P) values for the association test from the GWAS meta-analysis, and the x-axis the chromosomal positions (GRCh38/hg38). The SNP with the lowest P value in the locus is indicated by a purple diamond.
The remaining SNPs in the region are colored according to their r2 with the lead SNP (r2 calculated using 1000 Genomes European (EUR) reference). The dotted gray line represents the genome-wide significance threshold of 5 × 10−8.Table 2Association of discovered loci and independent SNPs with EFS (overall, 0–24 months and >24 months) in follicular lymphoma.Entire time period≤24 months>24 monthsCytobandNearest gene(s)SNP rsIDREFEAMAFCohortsNEventsHR (95% CI)P valueQPhetI2 %LOCONEventsHR (95% CI)P valueQPhetI2 %LOCONEventsHR (95% CI)P valueQPhetI2 %LOCO3q21.1KALRNrs1440146TC0.42MER2561341.46 (1.13–1.88)3.84E-030256611.86 (1.25–2.77)2.08E-030195731.21 (0.87–1.69)2.58E-010PRIMA3862091.34 (1.10–1.62)3.31E-031386871.56 (1.15–2.11)4.46E-0312921221.20 (0.94–1.54)1.41E-010FOLL05197941.27 (0.94–1.73)1.24E-010197521.42 (0.93–2.16)1.07E-010144421.13 (0.73–1.76)5.83E-010RELEVANCE203871.33 (0.95–1.88)9.94E-020203381.63 (0.99–2.69)5.35E-021162491.13 (0.72–1.77)5.93E-010POOLED10425241.35 (1.19–1.53)2.24E-060.500.9201/410422381.61 (1.33–1.95)9.90E-072.760.2527.662/47932861.18 (1.00–1.39)4.97E-022.760.2527.660/43q21.1KALRNrs333302TC0.42MER2531331.45 (1.12–1.87)4.69E-030253611.87 (1.26–2.77)1.96E-030192721.19 (0.85–1.66)3.08E-010PRIMA3832091.32 (1.09–1.61)4.65E-030383871.53 (1.13–2.08)5.92E-0312901221.20 (0.93–1.53)1.60E-010FOLL05196941.31 (0.96–1.78)9.02E-021196521.45 (0.95–2.21)8.79E-021143421.16 (0.75–1.81)4.98E-010RELEVANCE202871.33 (0.94–1.87)1.02E-010202381.63 (0.99–2.67)5.47E-021161491.13 (0.72–1.76)6.00E-010POOLED10345231.35 (1.19–1.53)3.07E-060.950.3801/410342381.61 (1.33–1.95)1.15E-062.480.2919.323/47862851.17 (1.00–1.39)5.72E-022.480.2919.320/43q27.3HRG/FETUBrs72625024TC0.09MER2531330.54 (0.38–0.77)5.56E-041253600.49 (0.31–0.77)1.78E-030193730.62 (0.36–1.07)8.77E-021PRIMA3802050.67 (0.50–0.89)5.92E-031380870.69 (0.45–1.07)9.82E-0202861180.65 (0.44–0.94)2.37E-021FOLL05192890.60 (0.37–0.95)2.78E-021192490.82 (0.43–1.59)5.65E-010142400.41 (0.22–0.77)5.83E-030RELEVANCE199870.71 (0.48–1.06)9.50E-021199380.82 (0.44–1.52)5.34E-010158490.64 (0.39–1.06)8.39E-020POOLED10245140.62 (0.52–0.74)9.37E-081.270.7404/410242340.67 (0.52–0.86)1.63E-033.980.1449.800/47792800.59 (0.46–0.74)9.88E-063.980.1449.802/48p23.2CSMD1rs6558733CT0.19MER2591351.55 (1.08–2.23)1.68E-021259622.09 (1.13–3.87)1.89E-020197731.29 (0.83–2.02)2.59E-010PRIMA3882081.42 (1.08–1.88)1.30E-020388881.21 (0.80–1.83)3.68E-0102931201.61 (1.11–2.33)1.27E-020FOLL05198951.65 (1.10–2.49)1.60E-020198521.44 (0.84–2.48)1.82E-010145431.94 (1.04–3.59)3.63E-020RELEVANCE205891.57 (1.01–2.44)4.67E-021205391.11 (0.62–1.96)7.31E-010163502.26 (1.14–4.46)1.89E-020POOLED10505271.50 (1.25–1.78)7.44E-060.420.9402/410502411.38 (1.07–1.79)1.44E-022.710.2626.130/47982861.60 (1.26–2.03)1.27E-042.710.2626.130/412q23.3WSCD2/CMKLR1rs12824955TC0.27MER2581350.68 (0.51–0.91)1.01E-020258620.51 (0.34–0.76)8.01E-041196730.94 (0.61–1.46)7.87E-010PRIMA3882080.73 (0.58–0.91)4.31E-030388870.59 (0.43–0.80)9.28E-0412941210.87 (0.65–1.19)3.90E-010FOLL05199950.85 (0.61–1.18)3.37E-010199520.72 (0.47–1.10)1.29E-010146431.06 (0.65–1.73)8.24E-010RELEVANCE205890.81 (0.58–1.13)2.15E-010205390.65 (0.40–1.06)8.55E-021163500.96 (0.61–1.49)8.42E-010POOLED10505270.75 (0.65–0.86)3.18E-051.290.7300/410502400.60 (0.49–0.72)1.46E-072.650.2724.633/47992870.93 (0.76–1.13)4.54E-012.650.2724.630/414q32.13TCL6rs114695031TC0.20MER2551340.87 (0.65–1.18)3.75E-010255610.64 (0.43–0.96)2.91E-021194731.22 (0.77–1.93)3.91E-010PRIMA3802040.75 (0.59–0.95)1.93E-020380870.51 (0.37–0.71)5.76E-0512861171.10 (0.76–1.57)6.17E-010FOLL05191920.72 (0.52–1.01)5.41E-020191500.58 (0.38–0.88)1.06E-021140420.98 (0.58–1.66)9.48E-010RELEVANCE203880.73 (0.52–1.02)6.30E-020203380.53 (0.33–0.85)8.66E-031162500.96 (0.59–1.56)8.75E-010POOLED10295180.79 (0.68–0.91)1.28E-030.960.8100/410292360.58 (0.47–0.70)1.93E-084.170.1252.084/47822821.11 (0.89–1.38)3.65E-014.170.1252.080/421q21.1PPIAP1rs118086521AG0.06MER2551330.33 (0.20–0.54)8.56E-060255610.32 (0.17–0.58)1.97E-040194720.36 (0.16–0.80)1.21E-020PRIMA3832060.86 (0.56–1.33)5.05E-010383860.57 (0.33–1.00)5.08E-0212901201.35 (0.68–2.66)3.92E-010FOLL05194910.91 (0.49–1.69)7.65E-010194490.58 (0.28–1.20)1.42E-010144421.99 (0.61–6.56)2.56E-010RELEVANCE202880.59 (0.32–1.11)1.01E-010202390.45 (0.22–0.91)2.61E-021160491.05 (0.32–3.47)9.34E-010POOLED10345180.64 (0.50–0.83)7.07E-0410.180.0270.50/410342350.46 (0.34–0.64)1.67E-062.440.3017.942/47882831.03 (0.67–1.60)8.87E-012.440.3017.940/4REF reference allele, EA effect (modeled) allele, MAF minor allele frequency, Q Cochran’s Q test, Phet Cochran’s Q test heterogeneity P value, I2 heterogeneity index.LOCO leave-one-cohort-out model; for each SNP, the total LOCO score was obtained by summing the scores obtained in the four LOCO rounds, ranging from 0/4 to 4/4 (see “Methods,” Fig. 2 and Supplementary Fig.
2).
SNPs associated with EFS24
Identification of SNPs specifically associated with events within the first 24 months (EFS24) and with events after 24 months was performed using Cox modeling [22] piecewise. In the analysis for EFS24 (QQ plots, Supplementary Fig. 6 and Manhattan plots, Fig.
3b, c), four regions achieved a P < 10−5 with one region at 14q32.13 (rs114695031, P = 1.93 × 10−8) predicting events within the first 24 months at GWS (Supplementary Table 4). rs114695031 at 14q32.13 was discovered and validated in all four LOCO rounds and yielded an HR = 0.58 (P = 1.93 × 10−8) in the pooled analysis (Fig. 3b; Table 2 and Supplementary Table 4). This SNP is in the T cell Leukemia/Lymphoma 6 (TCL6) gene (Fig.
4d). Two suggestive SNPs associated with EFS24 localized at 12q23.3 (rs12824955 HR = 0.60; P = 1.46 × 10−7) and at 3q21.1 (rs333302 HR = 1.61; P = 1.15 × 10−6) had consistent results across cohorts with validation in three of four LOCO rounds (Fig. 3b, Table 2 and Supplementary Table 4).
The SNP rs12824955 at 12q23.3 is an intergenic SNP near WSC domain containing 2 (WSCD2) and Chemokine-like receptor 1 (CMKLR1) genes (Fig. 4e) and SNP rs333302 at 3q21.1 is an intronic SNP in the Kalirin RhoGEF kinase (KALRN) gene (Fig. 4f).
Finally, the SNP rs118086521 was associated with EFS24 (HR = 0.46; P = 1.67 × 10−6) (Table 2), localized at 21q21.1 near Peptidylprolyl Isomerase A Pseudogene 1 (PPIAP1) (Fig. 4g), and was validated in 2 of 4 LOCO rounds (Fig. 3b, and Supplementary Table 4).
Associations with events after 24 months (conditioned on surviving 24 months) showed no loci at GWS; there were three loci at P < 10−5 also discovered in the EFS analysis on the entire time period (Supplementary Table 5). Of note, sensitivity analyses in the MER cohort for top SNPs showed no difference when EFS was defined from the date of diagnosis or initiation of treatment (data not shown).
Polygenic risk score for EFS
A PRS was calculated using the top 6 independent SNPs from Table 2; rs1440146 and rs333302 were in strong LD, so we used rs333302 in the PRS. The 6-year EFS rates were 68% (95% CI: 62–75), 68% (95% CI: 62–76), 61% (95% CI: 55–69), 54% (95% CI: 48–62), and 31% (95% CI: 25–38) across quintiles of the 6-SNP PRS (Fig. 5a).
Analyses by individual cohorts (PRIMA, FOLL05, MER, RELEVANCE) are shown in Supplementary Fig. 7. Evaluation of EFS according to PRS by FLIPI score (0–2 and 3–5) showed good risk discrimination (Fig.
5b, c and Supplementary Table 6). After adjustment for sex, FLIPI, type of IC, cohort, and use of RM (modeled as a time-dependent variable), the PRS as a continuous variable had an HR = 1.93 (P = 3.41 × 10−25). Using Q3 as the reference, Q1 had an HR = 0.76 while Q5 had an HR = 2.56 (Supplementary Table 7).
Of note, for these analyses, we used a “saturated” model using all information for the FLIPI score (0, 1, 2, 3, 4, 5) in PRIMA, MER, and RELEVANCE, and in two categories (0–2 and 3–5) for FOLL05 (only grouped data available). In addition, using Harrell’s concordance index [26], the model with FLIPI and PRS had a better predictive ability for EFS than the model with FLIPI alone (c-statistic corrected by cross-validation, 0.67 [0.65–0.70] versus 0.62 [0.59–0.64]) (Supplementary Table 8).Fig. 5Kaplan–Meier survival curves according to the 6-SNP polygenic risk score (PRS).Event-free survival for the entire time period for the pooled cohorts (PRIMA, RELEVANCE, FOLL05, MER) (a); event-free survival for the entire period for the pooled cohorts stratified by FLIPI 0–2 (b) and FLIPI 3–5 (c); overall survival for the entire time period for the pooled cohorts (d); and stratified by FLIPI 0–2 (e) and FLIPI 3–5 (f).
Top SNPs associated with OS and HT
The top SNPs associated with EFS were analyzed for their association with OS. With the limitation of the number of events related to OS (Table 1), results obtained with EFS were also observed for OS (Supplementary Table 9). Patients in the highest quintile of the 6-SNP PRS (Q5: 2.41, 4.58) had a 6-year OS rate of 80% as compared to ≥90% for all other quintiles (Fig.
5d). Furthermore, patients with a FLIPI 0–2 and 3–5 with a high PRS had a 6-year OS rate of 86% and 71%, respectively (Supplementary Table 6 and Fig. 5e, f).
As HT has a major impact on FL prognosis and could arise during the initial period of 2 years following the beginning of IC suggesting an occult HT at diagnosis [32], we conducted an exploratory analysis in MER and PRIMA of the top SNPs modeled as competing risk of progression, HT or death as first event (Supplementary Fig. 8). With the limitation of the small number of events HT events (27 in MER and 17 in PRIMA), we found that rs12824955 at 12q23.2 and rs118086521 at 21q21.1 were associated with both HT and progression without HT, while rs6558733 at 8p23.2 and rs72625024 at 3q27.3 and rs333302 and rs1440146 at 3q21.1 were associated with only progression without transformation (Supplementary Fig.
9).
Multi-SNP PRS for EFS in cohorts with other first-line treatments
We next tested the prognostic value of the six SNP PRS discovered in IC-treated patients in three independent cohorts with other treatments. These cohorts were FL patients treated by rituximab and bendamustine in the MER (N = 82), rituximab monotherapy in the MER (N = 102), and rituximab in combination with lenalidomide in RELEVANCE (N = 169), corresponding to the experimental arm in this trial. We found that the prognostic effect of PRS discovered for IC patients was not predictive in these three cohorts, with the limitation of small numbers of patients and events compared to the IC-treated cohort (Supplementary Fig.
10a–c).
Functional assessment of target SNPs using tumor transcriptomics
To identify functional consequences of the SNPs associated with EFS, we used RNAseq data generated from 222 FFPE diagnostic tumors and paired GWAS from the RELEVANCE cohort to correlate the genotypes of our top SNPs with expression of coding genes (Gencode v24 as reference) within 500 kb of each sentinel SNP (Supplementary Table 10). While overall no significant results were observed after correction for multiple testing, several nominal associations were observed. The effect allele of rs72625024 (near HRG and FETUB and associated with superior EFS) correlated with lower expression of FETUB and RP11-573D15.9 and higher expression of RP11-78H24.1 (P = 0.20 for each coding gene) (Supplementary Fig.
11a). RP11-573D15.9 and RP11-78H24.1 are affiliated with the long-non-coding (lnc) RNA class. The effect allele of rs114695031 (near TCL6 and associated with superior EFS24) correlated with higher expression of SCARNA13 (Small Cajal Body-Specific RNA 13) (P = 0.11) (Supplementary Fig.
11b), known to mediate proliferation, invasion, and migration in a hepatocellular carcinoma model [33].
Functional assessment of target SNPs using open chromatin region (OCR) profiles
We then used ATACseq data from 56 diagnostic tumors and paired GWAS of patients from the RELEVANCE cohort to correlate the genotypes of our top SNPs with OCR profiles. For these exploratory analyses, we analyzed index SNPs identified in the present GWAS and SNPs in high LD (r2 = 1 and D’ = 1 in the 1000 Genomes individuals of European ancestry) and analyzed if variation of ATACseq peaks was observed at each position. For rs72625024 at 3q27.3 (near HRG and Fetuin), we found that this SNP was localized in active chromatin sites (Supplementary Fig.
12a) but with no statistical variation among ATACseq peaks between the three genotypes of the index SNP or SNPs in LD (Supplementary Fig. 12b). No peak was observed for rs114695031 at 14q32.13 (near TCL6) (Supplementary Fig.
12c).
Functional analyses from an in silico bioinformatics database
HaploReg and RegulomeDB analyses using data from primary B-cells data provided some evidence that index SNP rs72625024 near FETUB may have regulatory functions based on transcriptional factor binding sites, DNAse hypersensitivity, epigenetic marks, and distal enhancer-like cis-regulatory elements (Supplementary Fig. 13a–c). Similar findings were observed for the top EFS24 SNP rs114695031 near TCL6 with enhancer histone marks and motifs changed, with multiple additional SNPs in LD linked with promoter histone marks and DNAse hypersensitivity (Supplementary Fig.
13d–f). Global bioinformatics data support functionality of rs72625024 near FETUB with a global ranking score 1f, a probability score of 0.71 using RegulomeDB (Supplementary Fig. 14a) and 8 using FORGEdb (Supplementary Fig.
14b). Global ranking score is 2b for rs114695031 near TCL6 with a probability score of 0.61 and 5 using RegulomeDB and FORGEdb, respectively. A summary of RegulomeDB and FORGEdb scores for all top SNPs is provided in Supplementary Table 11.
Discussion
This GWAS of prognosis in patients with FL treated with first-line IC identified loci that were associated with EFS independent of FLIPI, type of IC, and RM. These results were generated from patients in three randomized phase 3 trials and one prospective real-world cohort. All studies had high-quality pathology review, key clinical data, and relevant clinical outcomes.
Outcome data were mature with a minimal median of follow-up for RELEVANCE of 81 months and comparable EFS across the four cohorts (6-year EFS between 54 and 60%), allowing the ability to assess both early clinical failure (as assessed by EFS24) as well as longer-term outcomes.
The locus at 14q32.13 yielded GWS for the prediction of early events (EFS24) in FL patients. We also identified five loci that were suggestive of an association with EFS (3q27.3, 8p23.2) and EFS24 (3q21.1, 12q23.3, 21q21.1). Further, a PRS constructed with these six independent SNPs improved the ability to predict EFS over FLIPI in multivariable analyses adjusted for type of chemotherapy and accounting for use of RM.
We acknowledge that since the PRS was developed and tested in the same dataset, it is at inherent risk of being over-fit. These results support the contribution of inherited genetic background in the prognostic evaluation of FL patients and the value of using agnostic genome-wide approaches. Indeed, in our evaluation of the prognostic effects of 36 SNPs reported in previous studies (candidate SNPs, one GWAS performed before the rituximab era, one large analysis of 6679 tag SNPs) (Supplementary Table 12), only one SNP (rs6025446) reaching nominal statistical significance (P < 0.05) for EFS and three SNPs (rs6025446, rs4806857, rs1421094, rs1801131) reaching it for OS (Supplementary Tables 13 and 14).
In our analysis, it was notable that some loci (14q32.12, 12q23.3, and 21q21.1) predicted only early events (i.e., with EFS24), and some were associated with risk over the entire time period (3q27.3 and 8p23.2). While we did not design our study to specifically evaluate the type of EFS event (i.e., FL progression vs. HT), two of the top SNPs (rs12824955 at 12q23.2 and rs118086521 at 21q21.1) were associated with both HT and progression without HT, while the other SNPs were mainly associated with progression without HT.
These exploratory results were based on a small number of transformations and require independent confirmation.
Our results suggest several potential new biological pathways relevant to FL prognosis. For the 3q27.3 locus near HRG/FETUB (marked by rs72625024), Fetuin-B, which is secreted from the liver and is a member of the type 3 cystatin gene superfamily, has shown tumor suppressor properties in several solid cancers via an oncogenic role via PI3K/AKT pathway [34]. HRG is also synthesized by the liver and regulates numerous biologic processes such as immune complex clearance, cell adhesion, angiogenesis, coagulation, and, consequently, tumor progression [35].
One study found that HRG expression in 22% of FL [36]. In a mouse model, HRG exerts several tumor suppressor activities through modulation of tumor-associated macrophages (TAM), especially the switch from TAM M2 to TAM M1 polarity profiles [37]. The 14q32.13 locus near TCL6 (marked by rs114695031) is affiliated to a long-non-coding RNA with tumor suppressor properties in some solid tumor models through SRC/AKT and PI3/AKT pathways [38].
PI3K/AKT and mTOR pathways are identified as an important oncogenic signaling in FL with the development of inhibitors [39]. In the recent study of Gao et al., GEP found that the PI3K-AKT–mTOR pathway was enriched in POD24 FL patients [40].
From a statistical perspective, the LOCO approach incorporates the concept of statistical decision-making based on probability thresholds with LOCO-train and test-sets to provide complementary information about the level of statistical significance of SNPs and cross-validation. More classical designs (e.g., with discovery and validation cohorts and then meta-analysis or a single pooled analysis) could lead to lack of reproducibility, particularly if the final results are driven by the discovery dataset corresponding to a “winner’s curse” phenomenon [41]. For instance, our pooled GWAS analysis showed suggestive loci associated with EFS but there was no cross-validation in any of the LOCO rounds (Examples in Supplementary Fig.
2). Of course, a limitation of our study is that we did not have an external validation cohort but our study had >75% statistical power for each of the four LOCO rounds to detect true findings for EFS and EFS24 at the GWS alpha level in a post-hoc power analysis (Supplementary Table 15).
One notable finding is that the PRS based on the top SNPs from IC-treated cohorts was not predictive of EFS for FL patients treated with other therapeutic strategies (R-bendamustine, rituximab-lenalidomide, rituximab monotherapy). This type result was previously observed for other prognostic biomarkers in FL, such as m7-FLIPI [5, 42], which was not prognostic for patients treated with R-monotherapy [43] and R-bendamustine [44]. Of course, this could be due to using a PRS that was over-fit and has limited performance in external datasets.
In summary, our GWAS found that loci at 3q27.3 and 14q32.13 were associated with EFS and EFS24, respectively, after adjustment for other prognostic factors and type of IC. Our study suggested that some germline SNPs may be associated specifically with early failures, which could underlie biologic pathways implicated in early relapses [40]. Our results should be confirmed in independent cohorts, with parallel efforts to understand functional biology, especially at single-cell resolution [45, 46] and to determine the impact of therapeutic strategies outside of IC.
Supplementary information
Supplementary InformationSupplementary Tables