Enhancer RNA transcription pinpoints functional genetic variants linked to asthma

Abstract

Bidirectional enhancer RNA (eRNA) transcription is a widespread response to environmental signals and glucocorticoids. We investigated whether single nucleotide polymorphisms (SNPs) within dynamically regulated eRNA-transcribing regions contribute to genetic variation in asthma. Through applying multivariate regression modeling with permutation-based significance thresholding to a large clinical cohort, we identified novel associations between asthma and 35 SNPs located in eRNA-transcribing regions implicated in regulating cellular processes relevant to asthma, including rs258760 (mean allele frequency = 0.34, asthma odds ratio = 0.95; P = 5.04E-03). We show that rs258760 disrupts an active aryl hydrocarbon receptor (AHR) response element linked to transcriptional regulation of the glucocorticoid receptor gene by AHR ligands, which are commonly found in combusted air pollution. The role of rs258760 as a protective variant for asthma was independently validated using UK Biobank data. Our findings establish eRNA signatures as a tool for discovery of functional genetic variants and define a novel association between air pollution, glucocorticoid signaling and asthma.

Fig. 1. Pipeline for discovery of novel asthma-SNP associations based on proximity to μ and examples of μ-SNP colocalization.

A Schematic of pipeline and filtering approaches to identify high confidence μ-SNPs associated with asthma. B–C GRO-seq tracks from BEAS-2B cells treated as indicated for 30 min and visualized in the Integrative Genomics Viewer (IGV) genome browser based on counts per million mapped reads (vertical scales). Blue indicates reads annotated to the sense strand while red indicates reads annotated to the antisense strand. The transcription start site (TSS) and direction of transcription are marked by arrows at the top of each screenshot. Magnified regions show the locations of μ originally calculated by Tfit (green) and following manual refinement (gray) relative to the indicated SNP within dynamically regulated enhancers containing high (B) or low (C) confidence μ-SNPs.

Fig. 2. High confidence μ-SNP regions recapitulate dynamic eRNA regulation patterns in reporter assays.

A Mean (±SD) normalized luciferase activity of indicated enhancer reporter constructs or empty vector (EV) control in BEAS-2B cells treated as indicated for 8 hr (n = 6 technical replicates per treatment). pFKBP5 and pTNFAIP3 were included as canonical positive controls for dex and TNF responses. *P_adj < 0.05 vs same reporter+veh; one-way ANOVA corrected for multiple comparisons using Bonferroni’s method (pFKBP5: P_dex < 0.0001, P_TNF+dex = 0.0017; pTNFAIP3: P_dex = 0.0008, P_TNF < 0.0001, P_TNF+dex < 0.0001; pCOL8A1: P_TNF < 0.0001, P_TNF+dex < 0.0001; pCFLAR: P_TNF < 0.0001, P_TNF+dex < 0.0001; pRHOB: P_dex < 0.0001, P_TNF+dex < 0.0001; pKCMF1: P_dex < 0.0001, P_TNF < 0.0001, P_TNF+dex < 0.0001; pANKRD1: P_dex = 0.0087; pCEBPB: P_dex < 0.0001, P_TNF < 0.0001; pERRFI1: P_TNF < 0.0001, P_TNF+dex < 0.0001; pIER3: P_dex = 0.0002, P_TNF < 0.0001, P_TNF+dex < 0.0001; pOSR1: P_dex < 0.0001, P_TNF < 0.0001, P_TNF+dex < 0.0001; pNR3C1: P_dex < 0.0001, P_TNF < 0.0001, P_TNF+dex < 0.0001; pGCNT1: P_dex < 0.0001, P_TNF = 0.0003, P_TNF+dex < 0.0001.) Source data are provided as a Source Data file. B GRO-seq tracks, as described for Fig. 1, for select regions interrogated by enhancer reporter assay.

Fig. 3. rs149411423 disrupts a functional glucocorticoid response element that regulates induction of GCNT1.

A GRO-seq tracks, as described for Fig.1, for the GCNT1 locus. Scaling applies through panel (D). B Glucocorticoid receptor (GR) ChIP-seq tracks from BEAS-2B (gray; top) and primary human airway smooth muscle (HASM) cells (black; bottom) ±1 hr dex visualized in the IGV browser based on counts per million mapped reads (vertical scales). C Micro-C chromatin contacts in BEAS-2B cells with arcs connecting 1 kb interacting regions across genomic space; darker color indicates higher contact frequency; blue indicates contacts of interest. Arc height is proportional to distance. D Heat map of chromatin contacts in publicly available Micro-C data generated in human embryonic stem cells; arrows indicate contact between regions of interest. E Magnified view of GRO-seq tracks within rs149411423-μ-SNP region shown in panel (A). Scaling applies through Panel (F). F Aligned ATAC-seq tracks from primary human airway epithelial cells cultured at air-liquid interface (ALI; green), based on counts per million mapped reads (vertical scales), and magnified view of GR ChIP-seq peaks described in panel (B). G Consensus binding logo (MatBase) for GR and sequence of match identified by MatInspector within the GCNT1 enhancer harboring rs149411423. H Mean (±SD) normalized luciferase activity of parent or point mutation GCNT1 reporter constructs in cells treated ±dex for 8 hr (BEAS-2B; n = 6 technical replicates except for pFKBP5 + dex, n = 5) or 24 hr (primary HASM; n = 6 technical replicates except for HASM1: pGCNT1+dex, n = 5 and HASM2: pFKBP5+veh, n = 5). *P_adj < 0.05 vs same reporter+veh; unpaired two-sided t-tests corrected for multiple comparisons using Holm-Sidak method (BEAS-2B: all P’s<0.000001; HASM1: P = 0.0028 for pFKBP5, P = 0.00005 for pGCNT1; HASM2: P = 0.0109 for pFKBP5, P = 0.0027 for pGCNT1, P = 0.0047 for pGCNT1 C > T), ^$P_adj < 0.05 vs EV+veh; one-way ANOVA corrected for multiple comparisons using Bonferroni’s test (BEAS-2B: P < 0.0001 for pGCNT1 and pGCNT1 C > T; HASM1: P = 0.0003 for pGCNT1; HASM2: P < 0.0001 for pGCNT1), ^#P_adj < 0.05 vs parent construct+dex; one-way ANOVA corrected for multiple comparisons using Bonferroni’s test (BEAS-2B: P < 0.0001 for pGCNT1 C > T; HASM1: P < 0.0001 for pGCNT1 C > T; HASM2: P < 0.0001 for pGCNT1 C > T). Source data are provided as a Source Data file.

Fig. 4. Expression and function of NR3C1 are regulated by aryl hydrocarbon receptor (AHR) ligands through an AHR response element disrupted by rs258760.

A IGV screenshots of BEAS-2B GRO-seq tracks aligned with PRO-seq tracks (sense strand light blue; antisense light purple) from BEAS-2B cells ±30 min wood smoke particles (WSP) treatment. B ATAC-seq tracks from primary human airway epithelial cells cultured at ALI. C Consensus binding logo (MatBase) for AHR and sequence match identified within the region containing rs258760. D Mean (±SD) normalized luciferase activity of indicated reporters (6 technical replicates) in BEAS-2B cells treated ±TCDD for 8 hr. ^P_adj = 0.05 vs EV+veh and *P_adj < 0.05 vs same reporter+veh; unpaired two-sided t-tests corrected here and throughout figure for multiple comparisons using Holm-Sidak test (P = 0.0039 for EV; P = 0.0003 for pMT2A, P < 0.0001 for pNR3C1, P = 0.0085 for pNR3C1 C > T), ^#P_adj < 0.05 vs parent construct⁺TCDD; one-way ANOVA corrected for multiple comparisons using Bonferroni method (P < 0.0001 for pNR3C1 C > T). E qPCR analysis of dex-mediated gene regulation following 4 hr TCDD pre-treatment in BEAS-2B cells (4 technical replicates). Bars depict mean (±SD) C_T values on a log₂ scale relative to DMSO+veh-treated cells. *P_adj < 0.05 vs same gene + DMSO+dex (unpaired two-sided t-tests corrected P = 0.0003 for FKBP5; P = 0.0006 for TNFAIP3; P = 0.0002 for NR3C1). F Magnified view of PRO-seq tracks from (A). G BEAS-2B Micro-C contacts as described for Fig. 3C (blue indicates relevant contacts). H Micro-C interactions as described for Fig. 3D; arrows indicate relevant contacts. I Visualization of SNPs by GWAS p-values for associations with forced vital capacity (FVC) across the chr5:143,100,000-143,400,000 (NR3C1) topologically-associated domain in the Lung Disease Knowledge Portal; association strength indicated by the color scale. J Mean (±SD) normalized luciferase activity of indicated reporters in BEAS-2B cells treated ±WSP for 8 hr (6 technical replicates). ^P_adj < 0.05 vs EV+veh and *P_adj < 0.05 vs same reporter+veh (unpaired two-sided t-tests corrected P = 0.0005 for EV, P < 0.0001 for pMT2A, P < 0.0001 for pNR3C1, P = 0.0015 for pNR3C1 C > T, P = 0.0005 for prs258753, P = 0.0081 for prs258753 T > C, P = 0.0005 for prs864354, P = 0.0002 for prs864354 G > T), ^$P_adj < 0.05 vs EV+veh; one-way ANOVA corrected for multiple comparisons using Bonferonni test (all P’s<0.0001), ^#p < 0.05 vs parent construct+WSP; one-way ANOVA Bonferonni corrected (P < 0.0001 for pNR3C1 C > T). Source data are provided as a Source Data file.