High-throughput identification of regulatory elements and functional assays to uncover susceptibility genes for nasopharyngeal carcinoma

Abstract

Large-scale genetic association studies have identified multiple susceptibility loci for nasopharyngeal carcinoma (NPC), but the underlying biological mechanisms remain to be explored. To gain insights into the genetic etiology of NPC, we conducted a follow-up study encompassing 6,907 cases and 10,472 controls and identified two additional NPC susceptibility loci, 9q22.33 (rs1867277; OR = 0.74, 95% CI = 0.68-0.81, p = 3.08 × 10^-11) and 17q12 (rs226241; OR = 1.42, 95% CI = 1.26-1.60, p = 1.62 × 10^-8). The two additional loci, together with two previously reported genome-wide significant loci, 5p15.33 and 9p21.3, were investigated by high-throughput sequencing for chromatin accessibility, histone modification, and promoter capture Hi-C (PCHi-C) profiling. Using luciferase reporter assays and CRISPR interference (CRISPRi) to validate the functional profiling, we identified PHF2 at locus 9q22.33 as a susceptibility gene. PHF2 encodes a histone demethylase and acts as a tumor suppressor. The risk alleles of the functional SNPs reduced the expression of the target gene PHF2 by inhibiting the enhancer activity of its long-range (4.3 Mb) cis-regulatory element, which promoted proliferation of NPC cells. In addition, we identified CDKN2B-AS1 as a susceptibility gene at locus 9p21.3, and the NPC risk allele of the functional SNP rs2069418 promoted the expression of CDKN2B-AS1 by increasing its enhancer activity. The overexpression of CDKN2B-AS1 facilitated proliferation of NPC cells. In summary, we identified functional SNPs and NPC susceptibility genes, which provides additional explanations for the genetic association signals and helps to uncover the underlying genetic etiology of NPC development.

Keywords: cis-regulatory elements; functional profiling; genome-wide association study; nasopharyngeal carcinoma.

Graphical abstract

Figure 1

Replication of GWAS meta-analysis identifies additional NPC risk loci (A) Study design. Published NPC meta-GWAS by He et al. (B) Additional SNPs associated with NPC risk. The odds ratios (ORs), 95% confidence intervals (95% CIs) and the corresponding p values for rs1867277 (9q22.33) and rs226241(17q12) were calculated by fixed-effect meta-analysis.

Figure 2

Integrated functional annotation for the seven SNPs at 9q22.33 locus (A) Chromatin interaction peaks identified by PCHi-C in C666-1 cells are shown in a looping format. (B) A partial enlarged view of the 9q22.33 locus. The SNPs are shown as red vertical lines according to their chromosomal positions (x axis). The seven potential functional SNPs at 9q22.33 co-localize with ATAC-seq peaks, histone (H3K4me1, H3K4me3, and H3K27ac) ChIP-seq peaks, and chromatin interaction peaks (yellow shading). The histone ChIP-seq peaks are shown as log likelihood ratio (called by MACS2) after subtracting the input of each cell line.

Figure 3

Integrated functional annotation for the two SNPs at 9p21.3 locus The SNPs are shown as red vertical lines according to their chromosomal positions (x axis). The two potential functional SNPs, rs2518723 and rs2069418, at 9p21.3 co-localize with ATAC-seq peaks, histone (H3K4me1, H3K4me3, and H3K27ac) ChIP-seq peaks, and chromatin interaction peaks (yellow shading). Chromatin interaction peaks identified by PCHi-C in C666-1 cells are shown in a looping format. The histone ChIP-seq peaks are shown as log likelihood ratio (called by MACS2) after subtracting the input of each cell line.

Figure 4

Experimental validation of potential functional SNPs at 9q22.33 locus (A) Luciferase activity of the DNA segments in the CRE1. The flanking sequence of the SNPs in CRE1 at 9q22.33 were cloned into the pGL3-promoter luciferase reporter vector and transfected into C666-1 and HK1 EBV+ cells. Five SNPs within CRE1 were mutated from the protective alleles to the risk alleles. Luciferase activity was normalized by pGL3 promoter (n = 5). (B) Two different CRISPRi single-guide (sg) RNAs (sgC1-1 and sgC1-2) were used to target dCas9-KRAB to CRE1. The sgCon contains a non-targeting control sgRNA. The expression of genes was measured by qPCR. The data in (A) and (B) are presented as mean ± standard deviation (SD). Student’s t test was used for statistical analysis. (C) eQTL analyses demonstrated the correlation between the genotypes of tag SNP rs1867277 and the expression of PHF2 in the 83 NPC tissues. The data are presented as mean and 95% confidence interval. Student’s t test was used for statistical analysis. (D) Proliferation curves before (vector) or after overexpression of PHF2 (PHF2) in C666-1 and HK1 EBV+ cells. Two-way ANOVA was used for the statistical analysis of cell viability assays. (E) Colony formation assay for C666-1 and HK1 EBV+ cells with PHF2 overexpression and control cells (vector). Student’s t test was used for statistical analysis. (F) Proliferation curves before (si-NC) or after knockdown of PHF2 (si-PHF2-1 and si-PHF2-2) in C666-1 and HK1 EBV+ cells. Two-way ANOVA was used for the statistical analysis of cell viability assays. (G) Colony formation assay for C666-1 and HK1 EBV+ cells with PHF2 knockdown (si-PHF2-1 and si-PHF2-2) and control cells (si-NC). Student’s t test was used for statistical analysis. The data in (D)–(G) are presented as mean ± SD. ^∗∗∗ indicates p < 0.001, ^∗∗ indicates p < 0.01, and ^∗ indicates p < 0.05.

Figure 5

Experimental validation of potential functional SNPs at 9p21.3 locus (A) The flanking sequence of the SNP rs2069418 of CRE3 was cloned into the pGL3-promoter luciferase reporter vector. Plasmids with different genotypes of rs2069418 were generated by site-directed mutagenesis and transfected into C666-1 and HK1 EBV+ cell lines. Luciferase activity results were normalized by pGL3 promoter (n = 5). (B) Two different CRISPRi single-guide (sg) RNAs (sgC3-1 and sgC3-2) were used to target dCas9-KRAB to CRE3 in C666-1 and HK1 EBV+ cells. The sgCon contains a non-targeting control sgRNA. The expression of genes was measured by qPCR. The data in (A) and (B) are presented as mean ± standard deviation (SD). Student’s t test was used for statistical analysis. (C) eQTL analyses demonstrated the correlation between the genotypes of rs2069418 and the expression of CDKN2B-AS1 in the 83 NPC tissues. The data are presented as mean and 95% confidence interval. Student’s t test was used for statistical analysis. (D) Proliferation curves before (vector) or after overexpression of CDKN2B-AS1 (CDKN2B-AS1) in C666-1 and HK1 EBV+ cells. Two-way ANOVA was used for the statistical analysis of cell viability assays. (E) Colony formation assay for C666-1 and HK1 EBV+ cells with CDKN2B-AS1 overexpression and control cells (vector). Student’s t test was used for statistical analysis. (F) Proliferation curves for C666-1 and HK1 EBV+ cells transfected with smart silencer of CDKN2B-AS1 (ss-CDKN2B-AS1) compared with their control groups (ss-NC). Two-way ANOVA was used for the statistical analysis of cell viability assays. (G) Colony formation assay for C666-1 and HK1 EBV+ cells transfected with smart silencer of CDKN2B-AS1 (ss-CDKN2B-AS1) compared with their control groups (ss-NC). Student’s t test was used for statistical analysis. The data in (D)-(G) are presented as mean ± SD. ^∗∗∗ indicates p < 0.001, ^∗∗ indicates p < 0.01, and ^∗ indicates p < 0.05.