Genetic variation near CXCL12 is associated with susceptibility to HIV-related non-Hodgkin lymphoma.

Human immunodeficiency virus (HIV) infection is associated with an increased risk of non-Hodgkin lymphoma (NHL). Even in the era of suppressive antiretroviral treatment, HIV-infected individuals remain at higher risk of developing NHL compared to the general population. To identify potential genetic risk loci, we performed case-control genome-wide association studies and a meta-analysis across three cohorts of HIV+ patients of European ancestry, including a total of 278 cases and 1924 matched controls. We observed a significant association with NHL susceptibility in the C-X-C motif chemokine ligand 12 (CXCL12) region on chromosome 10. A fine mapping analysis identified rs7919208 as the most likely causal variant (P = 4.77e-11), with the G>A polymorphism creating a new transcription factor binding site for BATF and JUND. These results suggest a modulatory role of CXCL12 regulation in the increased susceptibility to NHL observed in the HIV-infected population.

Introduction

Human immunodeficiency virus (HIV) infection is associated with a markedly increased risk of several types of cancer compared to the general population.[1-3] This elevated cancer risk can be attributed partly to viralinduced immunodeficiency, frequent co-infections with oncogenic viruses (e.g., Epstein-Barr virus [EBV], hepatitis B and hepatitis C viruses, human herpesvirus 8 [HHV- 8] and papillomavirus), and increased prevalence of traditional risk factors such as smoking.[4,5] However, all of these risk factors may not entirely explain the excess cancer burden seen in the HIV-infected (HIV+) population. [5,6]

A previous study performed in the Swiss HIV Cohort Study (SHCS) identified two AIDS-defining cancers, Kaposi sarcoma and non-Hodgkin lymphoma (NHL) as the main types of cancer found among HIV+ patients (NHL representing 34% of all identified cancers).[4] The relative risk of developing NHL in HIV+ patients was highly elevated compared to the general population (period-standardized incidence ratio [SIR] = 76.4).[4] High HIV plasma viral load, absence of antiretroviral therapy (ART) as well as low CD4+ T-cell counts are known predictive factors for NHL.[7,8] The introduction of ART into clinical practice has led to improved overall survival and restoration of immunity by decreasing viral load and increasing CD4+ T-cell counts, and has led to a decreased risk of developing NHL. However, the risk remains substantially elevated compared to the general population (suboptimal immune response [SIR] 9.1[range, 8.3–10.1])[9] and NHL still represents 20% of all cancers in people living with HIV in the ART era.[10] NHL associated with HIV are predominantly aggressive B-cell lymphomas. Although they are heterogeneous, they share several pathogenic mechanisms involving chronic antigen stimulation, impaired immune response, cytokine deregulation and reactivation of the oncogenic viruses EBV and HHV-8.[11]

The emergence of genome-wide approaches in human genomics has led to the discovery of many associations between common genetic polymorphisms and susceptibility to several diseases including HIV infection and multiple types of cancer.[12,13] Recent genome-wide association studies (GWAS) of NHL have identified multiple susceptibility loci in the European population.[14-22] These variants are located in the genes LPXN[21], BTNL2[23], EXOC2, NCOA1[14], PVT1[14,22], CXCR5, ETS1, LPP, and BCL2[22] for various subtypes of NHL, as well as BCL6 in the Chinese population.[24] Strong associations with variation in human leukocyte antigen (HLA) genes have also been reported.[15,18,22] However, in the setting of HIV infection, no genome-wide analysis has been reported concerning the occurrence of NHL and the specific mechanisms driving their development remain largely unknown.

Here we report the results of the first genome-wide analysis of NHL susceptibility in individuals chronically infected with HIV. We combined three HIV cohort studies from France, Switzerland and the USA and searched for associations between >6 million single nucleotide polymorphisms (SNP) and a diagnosis of NHL. We identified a novel genetic locus near CXCL12 to be associated with the development of NHL among HIV+ individuals.

Methods

Ethics statement

The SHCS, the Primo ANRS and ANRS CO16 Lymphovir cohorts (ANRS) and the Multicenter AIDS Cohort Study (MACS) cohorts have been approved by the competent Ethics Committees /Institutional Review Boards of all participating institutions. A written informed consent, including consent for human genetic testing, was obtained from all study participants.

Study participants and contributing centers

Genotyping and phenotypic data were obtained from a total of 2,202 HIV+ individuals enrolled in the SHCS, ANRS and MACS cohorts (278 cases and 1,924 controls) (Table 1). For details on inclusion criteria and cohorts, refer to the Online Supplementary Appendix.

Table 1.

Summary of included samples and studies.

Quality control and imputation of genotyping data

The genotyping data from each cohort was filtered and imputed in a similar way, with each genotyping array processed separately to minimize potential batch effects. All variants were first flipped to the correct strand orientation with BCFTOOLS (v1.8) using the human genome build GRCh37 as reference. Variants were removed if they had a larger than 20% minor allele frequency (MAF) deviation from the 1,000 genomes phase 3 EUR reference panel or if they showed a larger than 10% MAF deviation between genotyping chips in the same cohort.

The quality control (QC) filtered genotypes were phased with EAGLE225 and missing genotypes were imputed using PBWT[26] with the Sanger Imputation Service,[27] taking the 1,000 genomes project phase 3 panel as reference. Only high-quality variants with an imputation score (INFO >0.8) were retained for further analyses.

Genome-wide association testing and meta-analysis

In order to search for associations between human genomic variation and the development of HIV-related NHL, we first performed separate GWAS within each cohort (SHCS, ANRS and MACS) prior to combining the results in a meta-analysis.

For each cohort separately, the imputed variants were filtered out using PLINK (v2.00a2LM)[28] based on missingness (>0.1), MAF (<0.02) and deviation from Hardy-Weinberg Equilibrium (PHWE<1e-6). Determination of population structure and calculation of principal components was done using EIGENSTRAT (v6.1.4)[29] and the HapMap3 reference panel.[30] All individuals not clustering with the European HapMap3 samples were excluded from further analyses. The samples were screened using KING (v2.1.3)[31] to ensure no duplicate or cryptic related samples were included. Single-marker case-control association analyses were performed using linear mixed models, with genetic relationship matrices calculated between pairs of individuals according to the leave-one-chromosome-out principle, as implemented in GCTA mlma-loco (v1.91.4beta).[32,33] Sex was included as a covariate, except in the MACS cohort, which only includes men.

The results of the three GWAS were combined across cohorts using a weighted Z-score-based meta-analysis in PLINK (v1.90b5.4), after exclusion of the variants that were not present in all three cohorts.

Other methods

The details of the cohorts and other methods used, i.e., fine mapping, prediction of causal variants, long-range chromatin interactions, transcriptomic effects, comparisons to GWAS in the general population and information on data sharing can be found in the Online Supplementary Appendix.

Results

Study participants and association testing

In order to identify human genetic determinants of HIV-associated NHL, we performed case-control GWAS in three groups of HIV+ patients of European ancestry (SHCS, ANRS and MACS). The characteristics of the study participants are presented in Table 1. In total, genotyping data were obtained for 278 cases (NHL+/HIV+) and 1,924 matched controls (NHL-/HIV+). With this sample size, we had 80% power to detect a common genetic variant (10% minor allele frequency) with a relative risk of 2.5, assuming an additive genetic model and using Bonferroni correction for multiple testing (Pthreshold=5e-8).[34]

After genome-wide imputation and quality control, 6.2 million common variants were tested for association with the development of NHL using linear mixed models including sex as a covariate. Results were combined across cohorts using a weighted Z-score-based metaanalysis (Figure 1A). The genomic inflation factor (lambda) was in all cases very close to 1 [1.00–1.01], indicating an absence of systematic inflation of the association results (Figure 1B; Online Supplementary Figure S2).

Figure 1.

Genome-wide association analysis. (A) Schematic of analysis pipeline. (B) Quantile-quantile plot of the observed -log10(P-value) (black dots, y-axis) vs. expected -log10(P-values) under the null hypothesis (red line) to check for any genomic inflation of the observed P-values. No genomic inflation is observed, with the genomic inflation factor λ=0.99. (C) Manhattan plot of all obtained P-values for each variant included in the meta-analysis. The genome-wide threshold (P=5e-8) for significance is marked by a dotted line. Only variants at the CXCL12 locus were found to be significant.

Association results

We observed significant associations with the development of HIV-related NHL at a single locus on chromosome 10, downstream of CXCL12 (Figure 1C). A total of seven SNP in this locus had P-values lower than the genome-wide significance threshold (P<5e-8), with rs7919208 displaying the strongest association (Table 2). This association was only detected in the SHCS and ANRS cohorts and not among MACS study participants (Online Supplementary Table S1).

Table 2.

Significant association with human immunodeficiency virus-related non-Hodgkin lymphoma.

Fine mapping of the CXCL12 locus

In order to identify the causal variant(s) among associated SNP and determine their potential functional effects, we used a multi-level fine mapping approach, combining the statistical fine mapping tool PAINTOR to obtain a 99% credible set and the deep learning framework DeepSEA to predict any effects on chromatin marks and transcription factor binding these variants may have.

Using PAINTOR, we identified a single variant, rs7919208, having a high posterior probability (=100%) of being causal among the 99% credible set based on the integration of the association results, linkage disequilibrium (LD) structure and enrichment of genomic features in this locus (Figure 2).

Figure 2.

Fine mapping of genomewide significant hits with PAINTOR. (A) The 99% credible set and posterior probabilities of being the causal variant. The genomic positions are listed on the x-axis. Bottom tracks represent DNAase and chromatin marks obtained from GM12878 cells as well as transcription factor binding site (TFBS) from the Roadmap Epigenomics Project and ENCODE in the region. (B) Locus plot of the associated variants, highlighting the LD relationship, based on the Swiss HIV Cohort Study cohort. The top variant rs7919208 is marked by a black diamond.

Consistent with the PAINTOR result, DeepSEA also identified rs7919208 as the sole variant, among the 99% credible set, predicted to have a functional impact by significantly increasing the probability of binding by the Bcell transcription factors BATF (log2 fold-change=3.27) and JUND (log2 fold-change=2.91) (Online Supplementary Table S2). Further analysis of the genomic sequence surrounding rs7919208 and the JASPAR transcription factor binding site (TFBS) motifs for BATF and JUND revealed that rs7919208 G->A polymorphism creates the TFBS motif required for the binding of these transcription factors (Figure 3A).

Figure 3.

Novel transcription factor binding site and long-range interactions. (A) Canonical motifs of BATF and JUND with the underlying genomic reference sequence and the nucleotide change caused by rs7919208. (B) Promoter capture Hi-C analysis in the GM12878 cell line of the region with the predicted causal variant and CXCL12. Variants and their level of association in the meta-analysis are marked in the inner grey circle. Genome-wide significant variants are colored green. Purple lines indicate significant interactions between promoter and other genomic regions. (C) Topologically associating domains (TAD) in the GM12878 cell line in the region of CXCL12. The yellow and blue boxes indicate the called TAD from the Hi-C contact map above. The plot is centered on rs7919208. DHS: DNase I hypersensitive sites.

Long-range chromatin interactions

In order to assess the potential functional links between the TFBS created in the presence of the minor allele of rs7919208 and the nearby genes, we performed an analysis of promoter capture Hi-C data and topologically associating domains (TAD). We used the well-characterized GM12878 lymphoblastoid cell line produced by EBV transformation of B lymphocytes collected from a female European donor as a model.

First, in order to examine the interaction potential of the rs7919208 region with nearby promoters, we analyzed available promoter capture Hi-C data obtained from the GM12878 cell line. This analysis revealed a significant interaction between the rs7919208 region and the CXCL12 promoter, suggesting a possible modulating impact of rs7919208 on the transcription of that gene (Figure 3B). Second, in order to further validate this observed genomic interaction, we analyzed available TAD calls from GM12878 cells,[35] using the 3D Genome Browser for visualization[36] (Figure 3C). We observed that rs7919208 is located within a large TAD together with CXCL12, signifying the interaction potential of the new TFBS at rs7919208 and CXCL12.

Transcriptomic effects of rs7919208

We did not observe any association between rs7919208 and mRNA expression levels of CXCL12 in peripheral blood or peripheral blood mononuclear cells (PBMC) from multiple publicly available datasets, including GTEx,[37] GEUVADIS[38] and the Milieu Intérieur Consortium[39] (Online Supplementary Figure S3). Of note, CXCL12 expression levels were very low in all datasets (Online Supplementary Figure S4A).

Using allele-specific expression analysis in the GTEx dataset, we observed a significant effect in individuals heterozygous of rs7919208, with increased allelic imbalance of CXCL12 in fibroblasts (false discovery rate adjusted P=0.0006, one-sided rank sum test), which was not observed in other tissues (Online Supplementary Figure 4B).

HIV infection causes many profound transcriptomic changes.[40] Thus, in order to examine the effect of rs7919208 on CXCL12 in the context of HIV infection, we extracted RNA from PBMC of 452 individuals in the SHCS with available genotyping data and sequenced them using the Bulk RNA Barcoding and sequencing (BRB-seq) approach.[41] However, the expression levels of CXCL12 were below the limit of detection for most individuals, preventing an expression quantitative trait loci (eQTL) analysis.

Multiple isoforms of CXCL12 exist, with variable degrees of expression and potency described in the context of HIV infection.[42] We observed a single significant correlation between rs7919208 and CXCL12 transcript usage, which was restricted to visceral adipose tissue. The presence of the rs7919208 minor allele was associated with higher relative expression of the longer and rarer transcript isoform ENST00000374429.6 (Online Supplementary Figure S5).

No replication of susceptibility loci found in the general population

In order to assess whether the genetic contribution to the risk of developing NHL is similar or distinct in the HIV+ population compared to the general population, we extracted the P-values of all variants found to be genomewide significant in previous GWAS performed in the general population[14,21-24,43] and compared them to our results. We did not replicate any of the previously published genome-wide associated variants, even at nominal significance level (P<0.05), despite sufficient statistical power for many of the variants, thus indicating that the genetic susceptibility of NHL is distinct between the HIV+ and the general population (Online Supplemental Table S3). In order to further examine this possibility, we tested whether the NHL/HIV+ associated variant rs7919208 is associated with an increased risk of NHL in the general population. We performed a series of case/control GWAS of four NHL subtypes (CLL, DLBCL, FL and MZL) as well as a combined GWAS with all NHL subtypes (Online Supplementary Table S4; Online Supplementary Figure S6) and assessed the association evidence at rs7919208. We found no association between rs7919208 and any of the subtypes in the general population, even at nominal significance.

Discussion

In this genome-wide analysis, including a total of 278 NHL HIV+ cases and 1,924 HIV+ controls from three independent cohorts, we identified a novel NHL susceptibility locus on chromosome 10 near the CXCL12 gene. The strong signal observed in the meta-analysis was driven by the associations detected in the SHCS and ANRS cohorts and there was no evidence of association in the MACS cohort. Notably, most NHL cases in the MACS cohort date back to the pre-ART era, while only NHL cases diagnosed after the year 2000 were included in the SHCS and ANRS analyses. Conceivably, NHL occurring in the early years of the HIV pandemic may have been primarily driven by severe immunosuppression, which could have obscured any influence of human genetic variation among the cases in the MACS sample. Precise phenotype definition is crucial in designing large-scale genetic studies since any environmental noise tends to decrease the likelihood of identifying potential genetic influences.

NHL is a relatively rare cancer even among HIV infected individuals, making it difficult to collect the large numbers of cases that would typically be included in contemporary genome-wide genetic studies. Indeed, a recent study from the Data Collection on Adverse events of Anti-HIV Drugs (D:A:D) group showed an NHL incidence rate of 1.17/1,000 person-years of follow-up over the past 15 years (392 new cases in >40,000 HIV+ individuals). [8] Still, we were able to obtain clinical and genetic data from a total of 278 patients with confirmed NHL diagnosis. By matching them with a larger number of controls from the same cohorts, we had enough power to identify associated variants of relatively large effects in the CXCL12 region.

Several groups have already suggested a potential role for CXCL12 variation in HIV-related NHL. A prospective study correlated increased CXCL12 expression with subsequent NHL development in HIV+ children but not in uninfected children.[44] The number of A alleles at the CXCL12-3’ variant (rs1801157) has also previously been associated with an increased risk of developing HIV-related NHL during an 11.7 year follow-up period.[45] Thus, our data further support the role of CXCL12 as a critical modulator of the individual risk of developing NHL in the HIV+ population.

The role of CXCL12 and its receptor chemokine receptor 4 (CXCR4) in cancer in the general population is well established, with the levels of CXCL12 and CXCR4 found to be increased in multiple types of cancer and to be associated with tumor progression.[46,47] Furthermore, in vivo inhibition of either CXCR4 or CXCL12 signaling is capable of disrupting early lymphoma development in severe combined immunodeficient (SCID) mice transfused with EBV+ PBMC.[48] These results and others have already led to the development and testing of several small molecules targeting either CXCL12 or CXCR4 to inhibit tumor progression. [46]

We could not identify any significant relationship between rs7919208 and the expression levels of CXCL12 in PBMC or EBV transformed lymphocytes. This could be due to the low expression levels of CXCL12 in most tissues, apart from stromal cells. Still, our analysis of allelespecific expression showed a significant allelic imbalance for heterozygous carriers of rs7919208 for CXCL12. This signal was only observed in fibroblasts, the GTEx tissue most closely resembling stromal cells. Furthermore, we identified a significant association between rs7919208 and CXCL12 transcript isoform usage in the visceral adipose tissue, which is known to also contain a minority of stromal cells.[49] Combined these results underscore the potential importance of these cells in the development of HIV-related NHL and the ability of rs7919208 to modify transcription.

The new BATF and JUND binding site created by rs7919208 could act as an induced or dynamic eQTL, specifically triggered by HIV infection. Such eQTL can be found in regions deprived of regulatory annotations, since these have been mostly examined in static cell types.[50] Supporting this hypothesis, HIV has been shown to induce overexpression of BATF.[51] This would explain why rs7919208 is only a risk factor for HIV+ individuals and not in the general population.

Previous analyses in the general population have discovered both shared and distinct associations for NHL subtypes. [14,21-24,43] However, similar analyses were not possible in our sample since NHL subtype information was not available for many of our cases. Furthermore, information on serostatus for relevant co-infections with EBV or other oncogenic viruses was not available and could therefore not be assessed. In particular, EBV has been largely associated with the development of NHL and other lymphomas and is considered a driver of a subset of NHL in the general population.[52] Variants in the HLA region have consistently been associated with all NHL subtypes in HIV uninfected populations regardless of EBV serostatus, although different HLA associations have been observed for each NHL subtype. We did not find any evidence of HLA associations in our analyses of HIV-related NHL. This might be due to a lack of power, due to our limited sample size in comparison to the NHL GWAS performed in the general population. However, this lack of replication of HLA variants and all other risk variants previously identified in the general population strongly suggests that distinct genes or pathways influence susceptibility to NHL in the HIV+ population compared to the general population.[53] This distinction may be due to the unique pathogenic mechanisms involved in HIV-associated NHL, such as cytokine deregulation, chronic antigen stimulation and impaired immune response, among others.[11]

In summary, we have identified variants significantly associated with the development of NHL in the HIV+ population. Fine mapping of the associated locus and subsequent analyses of TAD, promoter capture Hi-C data as well as deep-learning models of mutational effects on transcription factor binding, points to a causative model involving the gain of a BATF and JUND transcription binding site downstream of CXCL12 capable of physically interacting with the CXCL12 promoter. These results suggest an important modulating role of CXCL12 in the development of HIV-related NHL.