Coronary Artery Disease Risk Variant Dampens the Expression of CALCRL by Reducing HSF Binding to Shear Stress Responsive Enhancer in Endothelial Cells In Vitro

The RNA-seq experiments reported in this study are deposited in the GEO (Gene Expression Omnibus) database under the accession number: GSE222118.

TeloHAEC (ATCC [American Type Culture Collection]; CRL-4052), HepG2 (ATCC; HB-8065), RAW 264.7 (ATCC; TIB-71), preadipocyte 3T3-L1 (ATCC; CL-173), and MOVAS (mouse aortic smooth muscle cells; ATCC; CRL-2797) cell lines were used in the luciferase experiment. TeloHAECs were cultured using vascular cell basal medium (ATCC; PCS-100-030) supplemented with Endothelial Cell Growth Kit-vascular endothelial growth factor (ATCC; PCS-100-041) and 10% fetal bovine serum (FBS). HepG2 cells were cultured in DMEM (4.5-g/L glucose, 2-mmol/L l-glutamine; Lonza) supplemented with 10% FBS (GIBCO; Thermo Fisher Scientific). RAW 264.7 was cultured using DMEM and FBS (heat inactivated in the water bath at 56° for 30 minutes). 3T3-L1 was cultured using DMEM (Lonza BE12-614F) supplemented with 10% FBS and 2-mM l-glutamine. MOVAS was cultured using DMEM supplemented with geneticin. Primary human coronary artery endothelial cells were cultured in media MV2 (purchased from Promocell) and used between passages 4 and 6. A concentration of 100-U/mL penicillin and 100-μg/mL streptomycin were used in all cell lines used in this study to prevent bacterial contamination.

Gene transcript coordinates were identified from GRO-seq (global run-on sequencing) data using a custom de novo detection pipeline.^15,16 Enhancer RNA (eRNA) quantification was performed using Homer¹⁷ analyzeRepeats.pl, with parameters -noadj, -noCondensing, and -pc 3. Quantification was done from the opposite strand of coding gene transcription for intragenic enhancers and from both strands for intergenic enhancers. Gene end coordinates were expanded by 1500 bp to account for transcription at the end of transcripts and annotating intragenic gene enhancers. Counts for eRNA and gene expression were normalized using the voom¹⁸ function of limma R package. Only enhancers with counts per million >1 in at least 5 samples were included in the analysis.

The elastic net multivariable regression model was applied to individual genes using eRNA as predictive variables. This approach was aimed at identifying the enhancers that are most effective in predicting gene expression. For each gene, only enhancers that were in the same topologically associating domain and P<0.05 from the Spearman correlation test of significance were used for fitting the model. Caret¹⁹ and glmnet²⁰ R packages and the repeatedcv method with 5 repeats and folds set to 5 were used for tuning and cross-validating the model. This process involved running the model across a range of α values from 0 to 1 in increments of 0.05 and recording the mean cross-validated error at the lambda.min value. Here, α and λ are regularization parameters. α controls the balance between the L1 (Lasso) and L2 (ridge) regularization terms. A higher α emphasizes fewer coefficients, while a lower α leans toward ridge regularization. Enhancers for a given gene were ranked based on their model coefficient to determine their importance for predicting gene expression. Model coefficients refer to the weights assigned to each feature in the predictive model. Therefore, when the model was fitted to predict gene expression values using enhancer expression values, the best enhancers got higher coefficients and were ranked higher than enhancers with coefficients close to zero.

The CRISPR reagents were adapted from the Alt-R system from IDT (Integrated DNA Technologies). The guide RNAs were designed using the IDT tool to minimize off-targeting effects using 2 guides to create 372-bp deletion (Table S1). The positive control (crRNA [CRISPR RNA] targeting HPRT [hypoxanthine phosphoribosyltransferase]) from the Alt-R clustered regularly interspaced short palindromic repeats/clustered regularly interspaced short palindromic repeat–associated 9 Human Control Kit (IDT) was used. Delivery of the ribonucleoprotein complex into TeloHAEC was performed using a Neon Transfection System with a 1350-V pulse of 30 ms width. Forty-eight hours after transfection, RNA was purified using an RNeasy Mini Kit (Qiagen) and the cDNA was prepared with a RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific). The mRNA level of the CALCRL was measured by SYBR green chemistry quantitative polymerase chain reaction (PCR) using specific primers (Table S2) in a StepOne Real-Time PCR System (Thermo Fisher Scientific). Three independent experiments with 4 technical replicates were performed. Data (ΔCt values) were checked for normal distribution using the Shapiro-Wilk test and the F test to compare variance before performing statistical tests, and the unpaired t test (2-tailed) was used. P<0.05 was used to define a significant difference between the groups. To generate clones (ΔEnh) with an enhancer deletion, cells underwent electroporation and were dispensed into a 96-well plate using a single-cell dilution technique. These cells were subsequently cultivated until clones exhibited the desired enhancer deletion.

scRNA-seq data previously generated from atherosclerotic human coronary arteries²¹ were obtained from the NCBI (National Center for Biotechnology Information) GEO repository (accession number GSE131778). The data processing and cell type annotation were performed as described in the study by Örd et al²² using the Seurat package.²³ For plotting gene expression in individual cells, the gene counts were depth-normalized to 10 000 total counts per cell and log-transformed.

DNA fragments containing the region rs880890-rs840585 with haplotypes G-T and A-C were amplified by PCR from genomic DNA with Phusion polymerase and specific primers (Table S3). For the dual-luciferase reporter assay, amplified DNA (698 bp) was subcloned into the hSTARR-seq_ORI plasmid (Addgene, No. 99296) that was used as a backbone for the plasmid constructs. The enhancer-inserted plasmid was cotransfected into TeloHAEC, HepG2, RAW 264.7, preadipocyte 3T3-L1, and MOVAS cell lines with the control vector pGL4.75 (Promega), which encodes the luciferase gene hRluc (Renilla reniformis). The luciferase constructs prepared with the inserts were verified by sequencing. Luciferase plasmids with minimal promoter containing rs840588(A)_rs840587(C), rs840588(G)_rs840587(A), rs907463(T), and rs907463(G) enhancers were purchased from vector builder. Transfection was performed in a 6-well plate using Lipofectamine Stem Transfection Reagent (Thermo Fisher Scientific; STEM00008) according to the protocol, and the Dual-Luciferase Reporter Assay System (Promega; E1980) was used to detect luciferase activity. Three independent experiments with 4 technical replicates were performed. Data followed a normal distribution (Shapiro-Wilk test) and equal variance (F test). Hence, the unpaired Student t test (2-tailed) was used for comparing enhancer activity in regions carrying rs880890 A allele compared with G.

TeloHAEC and HepG2 cell lines were used for the experiment. Self-transcribing active regulatory region sequencing library and experiment have been described previously by Örd et al,²² Selvarajan et al,²⁴ and Toropainen et al.²⁵ To assess the transcriptional activity of haplotypes of the candidate regulatory regions, self-transcribing active regulatory region sequencing RNA read data were UMI (unique molecular identifier)-deduplicated, depth-normalized to library size, and normalized for haplotype abundance in the plasmid DNA (input) library.

To identify the potential transcription factors binding to the region with rs880890, PERFECTOS-APE²⁶ was used with TFBS (transcription factor binding sites) motif collection, P value cutoff of 0.0005, fold change cutoff of 2.0, and motifbreakR²⁷ package with default settings.

Because TeloHAECs are homozygous to the region of the SNP rs880890, the following chromatin immunoprecipitation assay (ChIP)-seq data related to endothelial cells (HUVEC [human umbilical vein endothelial cells]) cultured under shear stress (for 6 hours at 37 °C and with constant CO₂ perfusion at ±3 dyn/cm²) and static were downloaded from GSE116241. Sequencing data were downloaded from NCBI using the SRA toolkit and processed using the nf-core ChIP-seq pipeline²⁸ to derive BAM files. The BaalChIP²⁹ tool was used for allele-specific measurements of transcription factor binding from the ChIP-seq data. A Bayesian statistical approach was used by the tool to correct the effect of background allele frequency on the observed ChIP-seq read counts. ATAC-seq (transposase-accessible chromatin with sequencing) data related to human aortic endothelial cells (HAECs) that were heterozygous at rs880890 undergoing shear stress were downloaded from GSE112340. Reads were aligned to hg19 with Bowtie2.³⁰ SAMtools³¹ was used to convert sam to BAM files. BAM files were sorted and indexed using the SAMtools sort and index function. IGV (integrative genomics viewer)³² was used to visualize the SNP position and count the read counts in unidirectional flow (UF) and disturbed flow (DF).

A flow device consisting of a computerized stepper motor UMD-17 (Arcus Technology) and a 1° tapered stainless steel cone was used to generate the physiologically relevant shear stress pattern.³³ The flow device was placed in a 37 °C incubator with 5% CO₂. The flow patterns used for this experiment can be found in Table S4. HAECs were seeded in 6-well plates at 100% confluence and cultured in flow media with EGM-2 (endothelial cell growth medium-2) medium (Lonza) containing 4% dextran. HAECs were subjected to UF or DF for 24 hours³³ before being harvested for RNA isolation. HAECs were seeded in 6-well plates in EGM-2 medium at 90% confluence and were transfected with 50-nmol/L siRNA using RNAiMAX (Life Technologies) as described by the manufacturer. siRNA targeting HSF1 (SI03246348) and nontargeting controls (1027310) were both purchased from Qiagen. Media was changed the following day to flow media with EGM-2 medium containing 4% dextran. HAECs were subjected to UF for another 24 hours³³ before being harvested for RNA isolation. RNA was isolated from cells using GenElute Mammalian Total RNA Miniprep Kit (Sigma-Aldrich) and reverse transcribed using a High-Capacity cDNA Reverse Transcription Kit (Thermo Fisher). Quantitative mRNA expression was determined by reverse transcriptase-quantitative PCR using SYBR Green MasterMix (Roche, Indianapolis, IN). The mRNA level of the HSF1 knockdown under shear stress was measured by using specific primers (Table S5).

For the luciferase reporter assay, the following flow setup was used. Ibidi Pump System Quad (10906; ibidi GmbH) and Ibidi µ-Slide VI with ibiTreat surface and channel height of 0.4 mm, length of 17 mm, and width of 3.8 mm (80606; ibidi GmbH) were used. Channels were first coated with PureCol collagen solution (5005; Advanced BioMatrix). Next, TeloHAECs (0.8×10⁵) were seeded into the channel in the vascular cell basal medium. After cell attachment, channels were rinsed once with medium, and perfusion sets (red, 10962; ibidi GmbH) were installed in the pump system. A shear stress of 10 dyn/cm² was first established for 30 minutes to allow the cells to adapt for flow, followed by an increase in shear stress to 20 dyn/cm² for 45 to 46 hours. For DF, an oscillatory flow setup for the Ibidi Pump System was used with 20 dyn/cm² and a switching time of 2 s (ie, 0.25 Hz). The flow rate remains constant except during valve switching simulating oscillatory turbulence.³⁴ The fluid velocity using the μ-Slide VI 0.4 10 dynes was 7.89 mL/min and 20 dynes was 15.77 mL/min. The viscosity (0.0072 dyns/cm²) of the fluid is a typical value for a cell culture medium supplemented with FBS at 37 °C. The formula used to calculate shear stress is given as follows: τ=η176.1Φ where Φ is the flow rate (mL/min), τ is the shear stress (dyn/m²), and η is the dynamical viscosity (dyns/cm²). For luciferase assay, channels were first rinsed twice with PBS, dry trypsinized, and channels were rinsed with lysis buffer to collect the cells. Data followed a normal distribution (Shapiro-Wilk test) and equal variance (F test). Hence, the unpaired Student t test (2-tailed) was used for comparing enhancer activity between rs880890_A and rs880890_G.

For the siRNA treatment, TeloHAECs were reverse transfected with the Lipofectamine RNAiMAX transfection reagent (Thermo Fisher) according to the manufacturer’s instructions. The target-specific Silencer Select siRNAs for CALCRL (s19889 and s19890) and scrambled siRNA were obtained from Thermo Fisher. Cells were lysed 24 hours after transfection. Three independent experiments with 4 technical replicates were performed. RNA extraction, cDNA synthesis, quantitative PCR, and statistical analysis followed the protocol described above for clustered regularly interspaced short palindromic repeats/clustered regularly interspaced short palindromic repeat–associated 9 experiments. The knockdown efficiency was >95%.

For the ChIP-quantitative PCR, the flow setup was the same as used in HSF1 Knockdown and CALCRL Expression Under Athero-Relevant Flows section. The flow device was placed in a 37 °C incubator with 5% CO₂. TeloHAECs at 100% confluence, maintained in EGM-2 medium containing 4% dextran (Sigma-Aldrich) in 6-well plates, were subjected to athero-protective UF or athero-susceptible DF for 24 hours before cells were processed for ChIP analysis. TeloHAECs subjected to UF or DF for 24 hours were cross-linked, digested, and immunoprecipitated according to the manual instruction of Pierce Agarose ChIP Kit (Thermo Scientific). Briefly, cells subjected to flow were washed with warm PBS before cross-linking with 1% formaldehyde. After a 10-minute incubation, glycine was added to 125-mmol/L final concentration, and, after 5 minutes, the cells were pelleted down and the supernatant was removed. Cells were then resuspended in 1-mL ice-cold PBS and 10-µL protease inhibitors (Halt Cocktail; Pierce) and were pelleted at 3000g for 5 minutes before lysing for undergoing MNase digestion at 10 U/µL for 15 minutes. After recovery of digested chromatin, the solution was immunoprecipitated with an anti-HSF1 antibody (Cell Signaling Technology; No. 4356) at 1.1 µg/2×10⁶ cells and with normal rabbit IgG as control overnight. After binding to agarose beads, the immunoprecipitate was eluted before DNA cleanup, and then, purified DNA was detected by quantitative PCR performed on LightCycler 480 II (Roche) using SYBR Green I Master (Roche) with primers (Table S6). Statistical analyses were performed with the 2-tailed unpaired Student t test.

TeloHAECs were collected (5×10⁶) in PBS by scraping from culture flasks and washed twice with cold PBS. Resuspension of cells was done into 500-μL 1× Hypotonic Buffer (20-mmol/L Tris-HCl, pH 7.4, 10-mmol/L NaCl, and 3-mmol/L MgCl₂) with protease inhibitor cocktail by pipetting up and down several times, followed by incubation on ice for 15 minutes. Twenty-five microliter of detergent (10% NP40) was added and vortexed for 10 s at the highest setting. Centrifugation was performed of the homogenate for 10 minutes at 3000 rpm at 4 °C. Cytoplasmic fraction in the supernatant was removed. The nuclear pellet was resuspended in 50-μL complete cell extraction buffer (10-mmol/L Tris, pH 7.4, 2-mmol/L Na₃VO₄, 100-mmol/L NaCl, 1% triton X-100, 1-mmol/L EDTA, 10% glycerol, 1-mmol/L EGTA, 0.1% SDS, 1-mmol/L NaF [sodium fluoride], 0.5% deoxycholate, and 20-mmol/L Na₄P₂O₇) for 30 minutes on ice by vortexing at 10-minute intervals followed by sonication, followed by centrifugation for 30 minutes at 14.000g at 4 °C. Quantitation of protein concentration was performed using the Quant-iT Protein Quantitation Kit (Thermo Fisher Scientific).

Oligonucleotide probes (Table S7; 15-bp flanking SNP site for reference or alternate allele) with a biotin tag at the 5′ end of the sequence (Integrated DNA Technologies) were incubated with TeloHAEC nuclear extract and the working reagent from the LightShift Chemiluminescent Electromobility Shift Assay Kit (Thermo Fisher Scientific, catalog number: 20148, Rockford, IL). For competitor assays, an unlabeled probe of the same sequence was added to the reaction mixture at 100× excess. The reaction was incubated for 30 minutes at room temperature and then loaded on a 6% retardation gel (Invitrogen, catalog number: EC6365BOX, California). The contents of the gel were transferred to a nylon membrane and visualized by the UV transillumination Image Lab Software (Bio-Rad Laboratories).

HAEC donors (n=2; Major Resources Table) heterozygous for the SNP rs880890 were selected to study the allele-specific binding of the variant to transcription factor HSF1. For ChIP, the experimental setup was the same as used in ChIP-Quantitative PCR Under the Athero-Relevant Flows section. Quantitative PCR was conducted using a custom TaqMan SNP Genotyping Assay for rs880890 (Thermo Fisher). This assay uses allele-specific probes coupled to different fluorescent dyes in the same reaction mix for the quantitative detection of the alleles in a single sample. The rs880890 assay amplifies a 130-bp region around this SNP. We calculated haplotype-specific enrichment for both anti-HSF1 and IgG. Fold enrichment was determined by dividing the enrichment of the immunoprecipitated chromatin by the enrichment of the respective input. The experiment was performed with 3 biological replicates, and quantitative PCR was performed in triplicate. Statistical analyses were performed with the 2-tailed unpaired Student t test.

Five micrometer of total cytoplasmic protein were combined with gel loading buffer, heated to 95 °C for 5 minutes and separated on Any kD Mini-PROTEAN TGX Stain-Free Protein Gels (No. 4568123; Bio-Rad Laboratories). Treated groups were loaded next to their corresponding controls at each time point in each gel. The proteins from enhancer-deleted clone (ΔEnh) and si-CALCRL were transferred to 0.2-µm nitrocellulose membranes (No. 1704158; Bio-Rad Laboratories) using the Trans-Blot Turbo Transfer System (Bio-Rad Laboratories), blocked in 5% nonfat dry milk in TBS-Tween and incubated for 1.5 hours at room temperature and with the primary antibody Anti-CALCRL (1:500, No. PA5-50644; Thermo Fisher Scientific) overnight. After washing, the membranes were incubated with horseradish peroxidase–conjugated anti-rabbit or anti-mouse secondary antibodies (1:2000, No. Sc-2357; Santa Cruz Biotechnology, or 1:1000, No. 7076; Cell Signaling Technology, respectively). The signal was visualized using Pierce ECL Plus Western Blotting Substrate (No. 32134; Thermo Fisher Scientific) and ChemiDoc (Bio-Rad Laboratories) for imaging. The membranes were stripped after probing with CALCRL primary antibody using RESTORE Western blot stripping buffer (No. PIER21059; Thermo Fisher Scientific) for 45 minutes at 37 °C. The membrane was then reprobed using the GAPDH (1:5000, No. 39-8600; Thermo Fisher Scientific) primary antibody. Three biological replicates were performed. The ratio of CALCRL to GAPDH was calculated using imageJ³⁵ for each lane, and the values of these ratios were normalized to the control group.

For the CALCRL-silenced and enhancer-deleted samples, total RNA quality was assessed using the Agilent 2100 Bioanalyzer System. RNA library preparation was handled using QuantSeq 3′ mRNA-seq Library Prep Kit FWD (Lexogen according to the manufacturer’s instructions). The library was quantified using Qubit dsDNA HS Assay Kit (Q32854; Thermo Fisher Scientific), and its quality was checked with a Bioanalyzer using High Sensitivity DNA Kit (5067-4626; Agilent Technologies). Individual libraries were pooled in equimolar ratio (4 nmol/L for each) and sequenced on an Illumina NextSeq 500. The bcl2fastq2 v2.20 (Illumina) was used to demultiplex sequencing data and convert base call files into FASTQ files. Nf-core²⁸ RNA-seq pipeline was used to process FASTQ files and derive gene counts. Gene counts were normalized for effective library size, and differentially expressed genes were analyzed using DESeq2.³⁶ Differentially expressed genes were defined by P<0.05 and an absolute fold change of >2.

Human coronary artery endothelial cells (purchased from Promocell) were seeded on 0.1% gelatin-coated slides with a density of 2.5×10⁵ and cultured for a minimum of 3 days to ensure complete confluence, production, and reorganization of subcellular matrix and maturation of cell-cell junctions. Culture underflow was performed for 72 hours using a parallel plate flow apparatus as described^37–39 under DF (5 dyn/cm²) or UF (1.5 dyn/cm²) shear stress. Total RNA was extracted from human coronary artery endothelial cells using the RNeasy Mini Kit (Norgen) according to the manufacturer’s protocol. RNA (100-ng) was reverse transcribed into cDNA with random primer by reverse transcriptase (Qiagen). cDNA (relative 1.7-ng RNA) was amplified by standard quantitative PCR with Taq DNA polymerase (SensiFAST, SYBR Green and Low ROX kit; Bioline). The geometric mean of GAPDH, RPLP0, and PABC4 quantification was used to normalize the expression of each test gene. The paired Student t test was used for data that followed normal distribution using the Shapiro-Wilk test. P<0.05 was used to define a significant difference between the groups (n=7 individual donors).

The weighted gene correlation network analysis⁴⁰ was used for the identification of coexpressed genes from RNA-seq samples. Weighted gene correlation network analysis is a well-established tool for studying biological networks based on pairwise correlation between variables in high-dimensional RNA-seq data sets. ComBat-seq⁴¹ was used to remove the batch effect from the data. Sequencing data were processed as mentioned in the RNA-Sequencing and Data Analysis section. After rlog transformation, we identified the genes, and samples with excessive numbers of missing samples were identified using the function goodSamplesGenesMS and were discarded from the analysis. The samples were then clustered using the function hclust to identify outliers in the data set. We constructed a weighted gene network using soft-thresholding power β. The function pickSoftThreshold was used, which aids in choosing a proper soft-thresholding power. We chose a power of 10, which resulted in an approximate scale-free topology network with the scale-free fitting index R2 >0.9. For the remaining genes, the choice of the soft-thresholding power β was used to which coexpression similarity is raised to calculate adjacency with function adjacency. The matrix of correlations was converted to an adjacency matrix of connection strengths. To minimize the effects of noise and spurious associations, the adjacency was transformed into a topological overlap matrix using TOM (topological overlap matrix) and the corresponding dissimilarity was calculated. The function hclust was used to hierarchically cluster the genes. Branches of the dendrogram group together densely interconnected are highly coexpressed genes. Modules were then defined as sets of genes with high topological overlap. The module hub gene was defined as the gene in the module with the highest connectivity or based on high intramodular connectivity. Cytoscape⁴² app geneMANIA⁴³ was used to identify known and novel coexpressed genes.

To identify how much of the coexpressed genes are enriched in the differential expressed data set of CALCRL depletion, a gene enrichment analysis was performed using the hypergeometric t test.⁴⁴ The parameters for the hypergeometric t test are given as follows. The number of successes (k=49) was the coexpressed genes significantly differentially expressed under CALCRL depletion (P<0.05), the sample size (s=195) was the total coexpressed genes, and the number of successes in the population (M=1460) represented all the genes that were significantly differently expressed under CALCRL depletion (P<0.05; in both siRNA and enhancer deletion), while population size (N=11 820) was the total number of genes included in the analysis.

Nitrate Colorimetric Assay Kit (Cayman Chemical; 780001) was used to measure the total nitrate production in an enhancer-deleted clone (ΔEnh). Cells from ΔEnh and control cells (from a similar passage as the enhancer-deleted clones) were plated in Ibidi µ-Slide VI and underwent unidirectional, disturbed, and static conditions as mentioned in Dual-Luciferase Reporter Assays Under Athero-Relevant Flows section. Culture media was collected after 48 hours. Total nitrate/nitrite concentrations of ΔEnh clone and control samples under UF and DF were normalized to static conditions. The assay was performed according to the manufacturer’s protocol. The LDH (lactate dehydrogenase) Cytotoxicity Assay Kit (Cayman Chemical; 601170) was used to measure cell death in the treatment of ox-LDL (oxidized low-density lipoprotein; Thermo Fisher; L34357) and LDL (Thermo Fisher; L3486). ΔEnh and control cells were plated in a 96-well plate and subjected to 200-µg/mL ox-LDL and LDL treatment for 24 hours. The assay was performed according to the manufacturer’s protocol. Statistical analyses were performed with the unpaired Student t test (2-tailed).

Transfected cells were seeded at 5×10³ cells per well in 200 µL of complete media in E-plates (ACEA Biosciences, San Diego, CA) and grown for 48 hours while monitoring with xCELLigence DP system (ACEA Biosciences) that measures electrical impedance across interdigitated gold microelectrodes integrated on the bottom of tissue culture plates. The xCELLigence recorded cell index readings every 1 hour for 3 days. Three biological replicates were used to determine cell proliferation and represented the relative numbers of cells compared with control cells. A 2-way repeated measures ANOVA test with the Sidak multiple comparisons test was used to compare CALCRL siRNA knockout treatment to scrambled (negative) control and ΔEnh clone to control, with P_adj≤0.05 deemed significant.

Forty-eight–well plate was coated with Matrigel basement membrane matrix 150 μL/well (Corning, growth factor reduced). For CALCRL knockdown and ΔEnh clones, TeloHAECs (40 000 cells) were seeded onto the Matrigel 48 hours after siRNA transfection. Knockdown efficiencies of >95% were consistently attained and detected by quantitative PCR. Images were taken from 9 spots per replicate every 4 hours until 24 hours using Incucyte (S3 System). Each siRNA knockdown had 3 technical replicates in each experiment, and experiments were repeated 3×. Analysis was done with ImageJ³⁵ and associated macro tool Angiogenesis Analyzer⁴⁵ where the total tube length was selected as the parameter of evaluation. For time-course analysis, each time point represents the average tube length from all spots of all replicates with the same treatments at the same time point. Data were checked for normal distribution before performing statistical tests. The unpaired Student t test (2-tailed) was used for data because values followed a normal distribution and equal variance.

CALCRL is a G-protein–coupled receptor and a recently identified CAD locus.¹¹ Further analysis of the SNPs within the CALCRL locus using the Open Target Genetics portal⁴⁶ demonstrated that the SNPs are also associated with hypertension (Figure 1A). Because the majority of the variants were located within noncoding regions and overlapped enhancer elements marked by DNase hypersensitivity and H3K27ac, we next sought to identify the relative importance of individual enhancers in the regulation of CALCRL expression. Starting from the premise that eRNAs are coexpressed with their target genes,^48,49 we constructed coexpression networks based on regularized logistic regression. This was achieved by comparing coactivity patterns of eRNAs and coding gene expression across cell types based on public GRO-seq data on 336 samples representing 45 cell types.¹⁵ Our model reports a linear combination of coefficients derived from eRNA expression to describe how important each enhancer is as a predictor for gene expression. Complex heatmap⁴⁷ of 12 regulatory regions (Figure 1B and 1C) was included in our analysis, and enhancer-8 demonstrated the highest regression coefficient of 0.31 and was, thus, predicted to be central for the regulation of CALCRL expression (Table S8). The plot of GRO-seq enhancer-8 eRNA and CALCRL mRNA expression showed a positive correlation with a P value of 1.38×10⁻¹¹ (Figure S1A). Hi-C (chromosome conformation capture) data from aortic endothelial cells, TeloHAECs,⁵⁰ further confirmed the interaction between enhancer-8 and the promoter of CALCRL (Figure S1B). Finally, to confirm the predicted regulatory role of enhancer-8 in controlling endothelial CALCRL expression, the Alt-R clustered regularly interspaced short palindromic repeats/clustered regularly interspaced short palindromic repeat–associated 9 system was used to selectively delete a 372-bp genomic region in TeloHAECs using electroporation. PCR assays and Western blot detected around 50% genetic deletion in the regulatory region enclosing rs880890 and a reduction in CALCRL expression in TeloHAECs (Figure S2). The enhancer-deleted cells (ΔEnh) exhibited significantly lower gene expression (P=0.003) compared with control cells confirming the importance of the enhancer elements in regulating CALCRL expression (Figure 1D).

The enhancer-8 harbored 5 CAD-associated SNPs, and we next focused on characterizing rs880890 and rs840585 in the DNAse hypersensitivity region in more detail. (Figure 1B, zoomed out). The GTEX⁵¹ (Genotype-Tissue Expression) database suggested that CALCRL expression is ubiquitous throughout tissues (Figure S3) with the highest expression seen in the lung, the artery, and the adipose tissue. To identify the cell type contributing to the expression of CALCRL in the arterial tissue, we took advantage of single-cell data from mouse single-cell atlas (Tabula Muris⁵²). tSNE (t-distributed stochastic neighbor embedding) visualization of all tissues from Tabula Muris identified CALCRL expression predominantly in endothelial cells (Figure 2A). In addition, analysis of public scRNA-seq data from human atherosclerotic lesions²¹ demonstrated high expression of CALCRL in endothelial cells and, to a lesser degree, in smooth muscle cells (Figure 2B; Figure S4). In line with the cell type–specific gene expression, enhancer-8 was only found accessible in endothelial cells of the human atherosclerotic aorta based on our previously published single-nucleus ATAC-seq data²² (Figure 2C).

To study whether the cell type specificity could be driven by enhancer-8, we measured the enhancer activity of 698-bp region in endothelial cells (TeloHAEC), preadipocytes (3T3L1), smooth muscle cells (MOVAS), hepatocytes (HepG2), and macrophages (RAW 264.7). Enhancer-8 harbored 5 CAD SNPs of which 2 were in the DNase hypersensitive region. The SNPs were assessed through a luciferase assay. While rs840588, rs840587, and rs907463 did not show any difference in enhancer activity (Figure S5A), haplotype (ref: rs880890[A]_rs840585[C] and alt: rs880890[G]_rs840585[T]) showed significant difference compared with the empty plasmid. Low enhancer activity was detected in all cell types except TeloHAECs, which demonstrated 100-fold higher activity. Importantly, when we compared the luciferase signal of the A allele with the G allele, there was a 40% increase in the activity of A compared with G (Figure 2D). This cis-eQTL was confirmed in the GTEx database where the rs880890 A allele was significantly associated with increased expression of CALCRL compared with the G allele (Figure S5B). Also supporting our results, a trend of allele-specific deposition of histone mark H3K27ac and the binding of endothelial-specific transcription factor ERG was also detected in a cohort of 21 to 44 HAEC donors⁵³ (Figure S6A). This is in line with the predicted ERG motifs in enhancer-8 situated 238 bp upstream and 291 bp downstream around the HSF1 motif at rs880890 SNP.

To separate between the 2 SNPs in mediating the allele-specific effect, we cloned the naturally occurring 3 haplotype combinations into a reporter vector with haplotype 1 containing the reference alleles rs880890(A) and rs840585(C) while haplotype 2 carrying rs880890(G) allele and haplotype 3 rs840585(T) allele. Our analysis suggests that compared with haplotype 1, only haplotype 2 showed allele-specific enhancer activity, while hap3 showed no effect (Figure 2E), pinpointing rs880890 as the most probable causal SNP. Importantly, rs880890 also demonstrated significant evolutionary conservation (PhyloP [phylogenetic P values] −log P=0.57) compared with rs840585 that was predicted to be fast evolving (PhyloP −log P=−0.99).⁵⁴ To understand the allele-specific binding of transcription factors to the rs880890, we performed an electromobility shift assay using biotinylated 31-bp probes targeting either the reference or alternative allele in TeloHAECs. Competitor assays were performed by incubating the reaction with ×100 excess of unlabeled (no biotin) oligonucleotide complexes with identical sequences (Figure S7). The results demonstrated an increased protein binding to the A allele compared with the G allele of rs880890 (Figure 3A). The sequence-based motif analysis revealed that HSFs such as HSF1, HSF2, and HSF4 are likely to bind to the CALCRL-regulating enhancer at the rs880890 SNP region with a preference toward the A allele in binding affinity (Figure S8; Table S9). HSFs are known to activate heat shock proteins during stress,⁵⁵ and in endothelial cells, HSF1 increase corresponds to elevated levels of eNOS (endothelial NO synthase).⁵⁶ To further explore this, allele-specific ChIP-quantitative PCR was conducted to assess HSF1 binding in heterozygous HAEC donors, using a custom TaqMan Genotyping Assay for rs880890. Confirming motif predictions, HSF1 showed a 1.63-fold higher enrichment at the A allele compared with the G allele (donor 1: P=0.004; donor 2: P=0.007; Figure 3B).

As CALCRL is a key enhancer of the eNOS pathway implicated in fluid shear stress response and HSF1⁵⁷ has been shown to be activated by shear stress, we wanted to investigate the effect of rs880890 SNP on enhancer activity under such conditions. As expected, endothelial cells exhibited morphological changes under UF, DF, and static flow conditions (Figure S9). Luciferase assay confirmed significantly lower enhancer activity after HSF1 knockdown, while an even further decrease in the luciferase activity of enhancers was observed under DF compared with UF (Figure 3C and 3D). ATAC-seq analysis of HAECs subjected to UF mimicking hemodynamics in the human distal carotid artery or the DF waveform mimicking the hemodynamics in the human carotid demonstrated that the rs880890-harboring enhancer is more accessible in HAECs under UF compared with DF (Figure 3E). Moreover, analysis of the publicly available H3K27ac ChIP-seq data allowed us to interrogate the allelic bias of rs880890 under shear stress in endothelial cells. The results confirmed that the A allele is indeed biased toward higher enhancer activity binding compared with the G allele, and this effect is more prominent under shear stress in HAECs and HUVECs (Figure 3F; Table S10). These effects were concordant with the regulation of CALCRL expression, which was higher in HAECs subjected to UF and lower in DF (P<0.001; Figure 3G). HSF1 knockdown under UF (Figure 3H) resulted in a significant decrease in CALCRL expression. This is in line with the directionality of HSF1 expression under shear stress with lower HSF1 levels detected under DF compared with UF.⁵⁸ In contrast, ERG knockout resulted in no significant change in CALCRL expression under UF, whereas we did note a significant decrease under DF. This suggests a potential coregulatory relationship between ERG and HSF1 (Figure S6B). To further validate the anticipated binding of HSF1 to enhancer-8, we conducted a targeted HSF1 ChIP-PCR. Our findings revealed a marked reduction in HSF1 binding to the enhancer region in HAECs under DF conditions compared with UF (Figure 3I), thus confirming that shear stress alters the binding of HSF1 to the enhancer.

To understand the downstream effects of CALCRL and CALCRL-regulating enhancers, we performed RNA-seq on siRNA-mediated CALCRL knockdown and CRISPR-mediated rs880890-enhancer deletion. We then used the rlog-transformed matrix of our RNA-seq data to construct a coexpression network using 11 820 genes across 22 samples. This allowed the identification of 31 distinct gene coexpression modules in our RNA-seq data depicted in distinct colors (Figure S10). The CALCRL gene was found in the module blue along with 143 coexpressed genes. A dendrogram of consensus module eigengenes obtained on the consensus correlation was used to merge the similar modules blue and cyan with a threshold of 0.1 resulting in the identification of 195 genes coexpressed with CALCRL. Selected genes are illustrated on a geneMANIA network (Figure 4A) that demonstrates the relationships between CALCRL and coexpressed genes. Interestingly, we observed potentially new CALCRL coexpressed genes, such as APLN (apelin), ADAMTS18, and AKAP12.

Further analysis of the differentially expressed genes identified 390 and 990 genes significantly regulated (P_adj<0.05) upon gene silencing or enhancer deletion, respectively (Tables S11 and S12). Importantly, this included 25% of the 195 CALCRL coexpressed genes, demonstrating significant enrichment (hypergeometric t test P=6.5×10⁻⁷). Functional enrichment analysis demonstrated that similar pathways were regulated by both the gene and enhancer repression (Figure 4B). The most important pathways affected included the repression of APLN, EDN1 (endothelin-1), ADM, eNOS, and angiopoietin signaling pathways. APLN and cyclooxygenase (PTGS2 [prostaglandin-endoperoxide synthase]) genes themselves were significantly downregulated in both RNA-seq experiments, indicating they could be direct downstream targets of CALCRL (Figure 4C and 4D). Among the differentially expressed genes, PTGS1 and PTGS2 were also downregulated in DF compared with UF further linking these genes to shear stress (Figure S11A and S11B). Interestingly, we also observed a downregulation trend of EDN1. The RNA-seq results also indicated that depletion of CALCRL expression leads to repression of TGF-β signaling (THBS1, PTGS2, and TGFBR2) and angiogenesis regulatory genes (PDPK1 and ANGPTL4). In line with this, knockdown of CALCRL and enhancer deletion led to a significant decrease in cell proliferation (Figure 5A and 5B) and tube formation (Figure 5C–5E).

Finally, we performed an LDH cytotoxicity assay on ΔEnh clones to assess the effect of LDL (200 µg/mL) and ox-LDL (200 µg/mL) on cell viability. The results demonstrated that the ox-LDL treatment induces a significant increase (P=0.026) in cell growth in the ΔEnh clone (Figure 5F) compared with control cells. To explore the difference in the production of NO under flow, ΔEnh clones and control cells were exposed to UF and DF followed by detection of NO production by using the Griess assay. In our results, as shown in Figure 5G and Figure S11C, we observe that ΔEnh clones under DF showed a significant reduction in NO production compared with UF (P=0.035) further confirming the vital role of enhancer-8 in the regulation of CALCRL for NO production.

Figures S1–S11

Tables S1–S12

Major Resources Table