Genetic effects on chromatin accessibility uncover mechanisms of liver gene regulation and quantitative traits

Ethics statement and liver tissue samples

Noncancerous liver tissue samples from deceased, deidentified human donors were previously obtained (Innocenti et al. 2011), and use of these liver tissue samples in the current study was considered to be non-human subject research by the institutional review boards at St. Jude Children's Research Hospital (Memphis, TN) and the University of North Carolina (Chapel Hill, NC). Cause of death was recorded for all but eight participants. The most common causes of death included cerebrovascular accident, head trauma, and anoxia.

Nuclei extraction and ATAC-seq library preparation

Nuclei were extracted from frozen liver tissue samples as previously described (Currin et al. 2021). We prepared ATAC-seq libraries using the Omni-ATAC protocol (Corces et al. 2017) with 5 µL Tn5 transposase per reaction, cleaned the transposase reaction and final libraries with Zymo DNA Clean & Concentrator (D4029), visualized and quantified the libraries using a TapeStation, and sequenced 150 bp paired-end reads on either an Illumina HiSeq 4000 or HiSeq X. We initially generated three ATAC-seq libraries and aimed for 75 million high-quality aligned read pairs per individual. We generated additional libraries and/or resequenced existing libraries to increase depth for some individuals (Supplemental Tables 2, 3). After quality control, the number of libraries per individual ranged from one to six, with most individuals having either three or four libraries.

ATAC-seq read alignment and quality control

We trimmed sequencing adapters and aligned ATAC-seq reads to the GRCh38 build of the human genome (International Human Genome Sequencing Consortium 2001) and generated high-quality (mapq > 20), nonduplicated alignments as previously described (Perrin et al. 2021). We converted BAM files to BED files using BEDTools (Quinlan 2014) and called peaks on each ATAC-seq library separately using MACS2 (Zhang et al. 2008) with the parameters -q 0.05 –nomodel –shift -100 –extsize 200 –keep-dup all. We removed ATAC libraries with fewer than 10 million raw read pairs, TSS enrichment less than four as calculated by ATAQV (Orchard et al. 2020), and <10% of high-quality alignments overlapping peaks from downstream analyses. We also removed two libraries that were identified as outliers using PCA of normalized read counts overlapping a previously defined set of liver ATAC peaks (see Supplemental Methods; Currin et al. 2021).

Validation of sample identity

Genotyping for 224 liver donors was previously described (Innocenti et al. 2011). We obtained genotypes from the NCBI Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) under accession number GSE26105. We tested for sample swaps between ATAC-seq and genotype data using verifyBamID (Jun et al. 2012) using variants overlapping existing liver ATAC peaks (see Supplemental Methods; Currin et al. 2021). We considered the best-matching ATAC library to a set of genotypes to be a confident match if chipmix <0.02. We corrected sample swaps and removed ATAC libraries from further analysis if they did not confidently match to genotypes from any of the 224 individuals.

Genotype quality control

We used KING v2.2.5 (Manichaikul et al. 2010) with the parameters –unrelated –degree 2 to remove highly related individuals with up to second-degree relatedness. Using PLINK v1.9 (Purcell et al. 2007), we removed genotypes with missingness across samples >5% and excluded individuals that had genotype missingness >5%, heterozygosity F coefficient greater than four standard deviations from the mean, and/or inconsistency between reported sex and sex estimated from genotypes (see Supplemental Methods).

Selection of the final sample set and identification of consensus ATAC peaks

After ATAC and genotype quality control, 477 ATAC libraries from 138 individuals remained. For each individual, we merged the BAM files of the high-quality libraries and removed duplicates using Picard MarkDuplicates (https://github.com/broadinstitute/picard) because some libraries for the same donor represent resequencing of the same library preparation. We then called peaks on merged libraries per individual using the same procedure that we used for single libraries (see section “ATAC-seq read alignment and quality control”).

We merged peaks across individuals using BEDTools (Quinlan 2014) and retained merged peaks found in at least seven individuals (5% of 138). To verify the reported sex of the ATAC libraries, we counted the number of ATAC reads overlapping merged peaks using featureCounts (Liao et al. 2014), calculated the proportion of sex chromosome peak counts found on the X Chromosome (Chr X counts/[Chr X counts + Chr Y counts]), and classified a sample as female if >99% of counts on the sex chromosomes were found on the X Chromosome and as male otherwise.

Merging peaks across individuals can result in wide peak boundaries and concatenation of multiple nearby accessible chromatin regions into one peak. Therefore, we developed a procedure to refine the boundaries of the merged peaks found in at least seven individuals. We first calculated the number of individuals that had a peak overlapping each base of each merged peak, which we define as coverage, using BEDTools coverage (Quinlan 2014). We converted the output of BEDTools coverage into a matrix of coverage values with peak IDs as rows and peak bases as columns. For each peak, we identified all local minima in coverage; peaks that contain only one peak should have only two local minima, whereas peaks containing multiple subpeaks will have three or more local minima. To prevent classifying small fluctuations in coverage as subpeaks, we only retained adjacent pairs of minima if the distance between them was at least 200 bp and if the maximum coverage between them was both at least seven and greater than twice the maximum of the two minima. To define consensus peak boundaries, we identified the position between the two local minima with maximum coverage and extended in both directions until we reached the outermost positions that had half of the maximum coverage. This is similar to the full width at half maximum procedure used to identify boundaries of DNase peaks in ENCODE (Meuleman et al. 2020).

We determined which consensus peaks were shared with ATAC peaks from our pilot liver caQTLs study (Currin et al. 2021) and liver tissue ATAC peaks from ENCODE (The ENCODE Project Consortium et al. 2020) using BEDTools (see Supplemental Methods; Quinlan 2014).

Genotype imputation

We performed genotype imputation for 157 individuals that passed genotype quality control and that had ATAC data using the TOPMed imputation server (Das et al. 2016; Taliun et al. 2021), selecting Eagle v2.4 (Loh et al. 2016) for phasing and Minimac4 v1.0.2 (Das et al. 2016) for imputation (see Supplemental Methods). We filtered the imputation results to contain only the 138 individuals that passed ATAC quality control and biallelic single-nucleotide polymorphisms (SNPs) or indels with imputation r² > 0.3 and minor allele count ≥ 10 using BCFtools v1.13 (Danecek et al. 2021).

Ancestry classifications and genotype principal components

We generated a VCF file containing the 138 liver donors and the 1000 Genomes (1000G) individuals and that contained 99,206 common variants in approximate linkage equilibrium (see Supplemental Methods). We performed two rounds of genotype PCA using PLINK (Purcell et al. 2007): one using only the 138 liver donors to use for QTL covariates and one using both the liver donors and 1000G individuals to use for ancestry inference. We used two methods to infer ancestry: weighted k-nearest neighbors (WKNN), implemented in the kknn R (R Core Team 2015) package, and support vector machine (SVM), implemented in the e1071 R package. For both methods, we used the first two genotype PCs to calculate distances between individuals and predicted the ancestry of liver donors using models trained on 1000G individuals. For WKNN, we estimated the optimal value of k to be four using the train.kknn function. Both methods agreed that 118 individuals were similar to European populations, 17 individuals were similar to African populations, and one individual was similar to admixed American populations; we considered the remaining two individuals for which the methods disagreed as similar to admixed populations.

caQTL mapping

We removed ATAC mapping bias using the WASP mapping pipeline (van de Geijn et al. 2015), quantified peak accessibility using featureCounts (Liao et al. 2014), adjusted for GC bias using EDASeq (Risso et al. 2011), and adjusted for library size and performed variance stabilization using DESeq2 (Love et al. 2014). We used ATAC PCs to control for unmeasured sources of variation. We mapped caQTLs using inverse-normal-transformed, variance-stabilized peak counts and variants within 1 kb of peak centers using FastQTL (Ongen et al. 2016) with the parameters –permute 1000 –normal –window 1000 –seed 123456789. We included sex, two genotype PCs, and 25 ATAC PCs as covariates. We included two genotype PCs because the percentage of variance explained leveled off after the first two PCs. We chose a window of 1 kb from peak centers to focus on variants close to and within the associated peaks. We chose 25 ATAC PCs because it resulted in the first local maximum of the number of significant caPeaks. To identify ATAC peaks with a significant caQTL (caPeaks), we estimated the false-discovery rate (FDR) of the FastQTL beta-approximated P-values for the most strongly associated variant per peak using the Benjamini–Hochberg (BH) procedure (Benjamini and Hochberg 1995) and classified caPeaks with FDR < 5% as significant. To test for long-range caQTL effects, we also ran FastQTL using variants within 1 Mb of peak centers with the parameter –window 1000000, keeping all other parameters identical to the 1 kb analysis. We chose a window of 1 Mb from peak centers to capture interactions between peaks at a locus but still within a typical distance for cis regulation.

We determined how many caQTLs from our previous study (Currin et al. 2021) were replicated in the current study (see Supplemental Methods). To compare the percentages of caPeaks and non-caPeaks (FastQTL beta-adjusted P-value > 0.5) overlapping different types of regulatory elements, we downloaded chromatin states lifted over to GRCh38 from the Roadmap Epigenomics Project (Roadmap Epigenomics Consortium et al. 2015; https://egg2.wustl.edu/roadmap/data/byFileType/chromhmmSegmentations/ChmmModels/core_K27ac/jointModel/final/all_hg38lift.dense.browserFiles.tgz) and assigned ATAC peaks to liver tissue (epigenome ID E066) chromatin states as previously described (Currin et al. 2021).

To identify caQTL that may be specific to non-European populations, we selected caQTL for which the lead variant was not polymorphic (MAF = 0) in donors similar to European populations. We also compared population-level lead variant MAF in our cohort to MAF in TOPMed (Taliun et al. 2021) populations present in the TOP-LD (Huang et al. 2022) resource.

Downsampling of caQTL sequencing depth and sample size

For sequencing depth, we downsampled the WASP-filtered BAM files to 20%, 40%, 60%, and 80% of the original depth using SAMtools (Danecek et al. 2021) with the -s parameter. We used five technical replicates for each downsampled percentage (random number seeds ranging from one to five). For each replicate, we quantified accessibility of the consensus peaks identified from the full-depth data and performed caQTL mapping using the same pipeline used for the full-depth data, including the caQTL covariates from the full-depth data.

For sample size, we generated five replicate random subsets of the full sample set corresponding to 20% (n = 27), 40% (n = 55), 60% (n = 82), and 80% (n = 110) of total sample size using R (random number seeds ranging from one to five) (R Core Team 2015). For each replicate, we generated a subset of the full-sample, variance-stabilized, peak-by-sample matrix containing the downsampled sample list, recalculated ATAC PCs, and performed caQTL testing using the peaks and variants used for the full sample list and including sex, two genotype PCs, and a variable number of ATAC PCs as covariates. We tested ATAC PCs from 0% to 25% of the downsampled sample size and selected the number of PCs that resulted in the first maximum of the number of significant caPeaks.

Identification of TF motifs disrupted by caQTL variants

To identify caQTL variants that may disrupt TF binding, we searched for cases in which the strength of a TF motif match changed between the alleles of a caQTL variant. Although ChIP-seq is a more direct measure of TF binding than motifs, measuring allelic differences in binding using ChIP-seq is limited by the number of TFs with available ChIP-seq data, and allelic imbalance in ChIP-seq reads can only be measured at heterozygous sites with sufficient read depth in the respective data sets. Consequently, we used TF motifs to more comprehensively survey allelic effects of caQTL variants on TF binding. We tested if caQTL variants were more likely to disrupt specific motifs among a set of 516 motifs from Cis-BP v1.02 (Weirauch et al. 2014) using a previously described protocol (Currin et al. 2021) with a few changes. To ensure that the set of background non-caQTL variants was larger than the set of caQTL variants, we did not use propensity score matching when selecting non-caPeaks (FastQTL beta-adjusted P > 0.5). However, we still adjusted for peak GC content and distance to the closest protein coding TSS in the logistic regression model that tested for motif disruption associated with caQTL status. We did not include peak width because it was not correlated with caQTL status. We tested for association of caQTL status with motif disruption for the 280 motifs that remained after selecting motifs disrupted by at least 20 caQTL variants and selecting the single motif per TF that had the most disruptions. We considered a motif to be significantly associated with caQTL status if the logistic regression P-value was less than 1.8 × 10⁻⁴ (0.05/280 motifs).

Identification of coordinately regulated caPeaks

We conservatively grouped caPeaks from the 1 Mb caQTL analysis whose lead variants were identical or in very high LD with each other (r² > 0.9, TOPMed Europeans) (Taliun et al. 2021; Huang et al. 2022). We defined distance between variants as the complement of LD values (1 – r²), performed single linkage hierarchical clustering using the hclust function in R (R Core Team 2015), and identified clusters as groups with a maximum distance of 0.1 (1–0.9).

Identification of driver peaks at coordinated caQTL signals

We used two complementary methods to predict driver peaks at coordinated caQTLs: proxy variant overlap and the PHM (Kumasaka et al. 2019), which tests for causal interactions between pairs of peaks. For the proxy variant overlap approach, we identified all proxy variants of all clumped lead variants of a caQTL signal (LD r² ≥ 0.8) and overlapped these proxies with all caPeaks from the 1 Mb analysis. We classified a coordinated caQTL as having a single driver peak if only one caPeak within the coordinated set overlapped proxies of the signal and no caPeaks within or outside the coordinated set overlapped proxies of the signal. We ran PHM to test for causal interactions between pairs of caPeaks from the 1 Mb analysis. We edited the bayeslm1.sh and bayeslm2.sh scripts provided with PHM to increase the –window-size parameter from 500,000 to 1,000,000 to test caPeak pairs up to 1 Mb apart. We classified a pair of peaks as having a causal interaction if the posterior probability that either peak regulates the other peak was greater than 0.5. We determined the number and identity of caPeaks regulating (upstream) and regulated by (downstream) each caPeak. To predict driver peaks at coordinated caQTLs, we searched for causal interactions between all caPeaks in the coordinated set, as well as any caPeaks outside the set that overlapped proxies of the coordinated caQTL signal. We considered a coordinated caQTL to have a single driver peak if only one of the tested peaks had at least one downstream peak within the coordinated set and no upstream peaks inside or outside of the coordinated set. If the proxy variant overlap and PHM methods disagreed on the driver peak for a coordinated caQTL, we classified the proxy variant overlap peak as the driver peak.

We tested for significant overlap between driver peaks identified by the two methods using a two-sided Fisher's exact test. The contingency table for the Fisher's exact test consisted of the number of coordinated signals with driver peaks predicted by both methods, proxy overlap only, PHM only, and neither method. To ensure that each coordinated signal was represented once, we assigned the driver peak to proxy overlap for the 46 coordinated signals in which PHM and proxy overlap disagreed on the driver peak prediction. We confirmed that the overlap between the two methods was still significant when retaining the PHM predictions at these 46 discordant signals (OR = 2.7, P = 4.6 × 10⁻²⁶).

Enrichment of TF ChIP-seq peaks and motifs

We downloaded a metadata file of liver TF ChIP-seq conservative IDR thresholded peaks from ENCODE (The ENCODE Project Consortium et al. 2020) and removed files that had any notes in the “Audit NOT_COMPLIANT” column. The resulting data set consisted of files for 17 TFs (Ramaker et al. 2017), 15 of which were mapped in two donors (ENCDO060AAA and ENCDO882mMZ); we merged ChIP-seq peak summits across multiple donors for the same TF. For driver peaks, we used the 1538 peaks from the combined set of PHM and proxy overlap predictions. For response peaks, we used the 3034 nondriver peaks from the 1538 coordinated sets with a predicted driver peak. We extracted the central 200 bases of driver and response peaks to mitigate enrichment biases owing to differences in peak width distributions between the driver and response peaks. We used BEDTools intersect (Quinlan 2014) to count the number of driver and response peaks overlapping the single base pair summits of TF ChIP-seq peaks separately per TF. We determined if driver peaks were significantly more likely to overlap ChIP-seq summits of each TF using a Fisher's exact test and classified enrichments with OR > 1 and P < 2.9 × 10⁻³ (0.05/17 TFs) as significant.

We tested for enrichment of 286 TF motifs in driver peaks relative to response peaks. We generated the set of 286 motifs by selecting a single motif per TF based on highest information content from the set of 516 motifs used in the caQTL motif disruption analysis. We scanned the central 200 bp of driver and response peaks for motif occurrences using FIMO (Grant et al. 2011) with the default P-value threshold of 1 × 10⁻⁴. We computed the number of unique driver and response peaks that overlapped each motif using BEDTools (Quinlan 2014), tested for enrichment using a Fisher's exact test, and classified motifs as significantly enriched in the driver peaks (OR > 1) or response peaks (OR < 1) if P < 1.8 × 10⁻⁴ (0.05/283 motifs with at least one occurrence in driver or response peaks).

Genome browser views of caQTL signals

We visualized the genome context of caQTL signals using Plotgardener v1.2.10 (Kramer et al. 2022). To make tracks of ATAC signal, we generated bedGraph files normalized for sample read depth using BEDTools genomecov (Quinlan 2014), converted the bedGraph files to bigWig files using UCSCTools (Hinrichs et al. 2006), imported the bigWig files into R with the Plotgardener readBigwig function (Kramer et al. 2022), and calculated the mean of signal across individuals for each genotype of the caQTL lead variant. We incorporated GENCODE v41 (Frankish et al. 2019) protein-coding and lncRNA genes into Plotgardener (see Supplemental Methods).

Predicting target genes of caPeaks

Identification of caPeak target gene categories

We used the Human Protein Atlas (Uhlén et al. 2015) to determine the tissue specificity of protein coding caPeak target genes (https://www.proteinatlas.org/download/proteinatlas.tsv.zip). We considered a gene to be enriched in liver tissue if the Protein Atlas considered the gene to be “tissue enriched,” “tissue enhanced,” or “group enriched” and if liver was one of the tissues in the group. We determined if caPeak target genes were more likely to be enriched in liver tissue relative to genes detected in the Protein Atlas that were not caPeak target genes using logistic regression. To identify target genes that are TFs, we selected TFs that had either a directly determined or inferred binding motif from Cis-BP v2.00 (Weirauch et al. 2014) and restricted to autosomal genes from GENCODE v41 (Frankish et al. 2019) by matching on Ensembl ID. We determined if caPeaks were more or less likely to be linked to TFs compared with non-TF genes by each of the four peak–gene links separately using logistic regression. For eQTL, we also ran an additional model providing the distance to the nearest TSS for the same gene as a covariate. Because not all genes linked to caPeaks were tested for liver eQTLs in GTEx, we reran the logistic regression models only using those genes. We also ran logistic regression models using gene expression level as a covariate. To compare gene expression levels between TF and non-TF genes, we used liver gene counts from GTEx (https://storage.googleapis.com/adult-gtex/bulk-gex/v8/rna-seq/counts-by-tissue/gene_reads_2017-06-05_v8_liver.gct.gz) and calculated counts per million (CPM) using edgeR (Robinson et al. 2010) with the trimmed mean of the m-values method to calculate library sizes. We defined the expression level of a gene as the median CPM across individuals. To identify caPeak target genes involved in drug response, we obtained a previously compiled list of genes involved in drug response (Etheridge et al. 2020) and restricted to autosomal genes present in GENCODE (Frankish et al. 2019) v41 with biotype of protein coding or lncRNA by matching on gene symbol.

Enrichment of coordinated peaks and eQTLs

We tested if caQTL signals with coordinated peaks were more likely to colocalize with GTEx liver eQTLs compared with caQTL signals for noncoordinated peaks using logistic regression. We tested additional models including promoter peak status or distance to nearest TSS as a covariate. For coordinated peaks, we considered the set to have a promoter peak if at least one peak was TSS proximal. We calculated TSS distance using the peak in the coordinated set closest to the nearest TSS. For most models, we considered a coordinated peak set to colocalize with an eQTL if the caQTL for at least one of the component peaks colocalized with an eQTL. As this approach could potentially overcount coordinated peak sets, we also tested models using the 931 driver peaks predicted from the proxy overlap method to test for eQTL colocalization, one peak per coordinated set. As the background peaks, we used noncoordinated caQTL peaks for which caQTL proxy variants (LD r² ≥ 0.8) overlapped the affected caPeak and no other caPeak. We also tested driver peak models adjusting for the absolute value of caQTL effect size.

Stratified LD score regression

We used stratified LD score regression (Finucane et al. 2015) implemented in the LDSC software to test if liver caPeaks were enriched for heritability of specific traits. We downloaded GWAS summary statistics for 703 UKBB traits with a heritability significance of z ≥ 4 generated by the Benjamin Neale laboratory (https://nealelab.github.io/UKBB_ldsc/downloads.html). We converted the coordinates of caPeaks from GRCh38 to GRCh37 using liftOver (Hinrichs et al. 2006) and ran LDSC with the GRCh37 reference files provided with LDSC, which were built from 1000G phase 3 European genotypes. We restricted analyses to HapMap SNPs as previously described (Finucane et al. 2015). For each trait, we ran a model using liver caPeaks and the baseline v1.2 model. We considered two measures of significance: the P-value of the coefficient, which measures the contribution to heritability specific to liver caPeaks after removing contributions from all data sets in the baseline, and the P-value of the enrichment FC (Enrichment_p), which includes contributions to heritability from liver caPeaks that are shared with the baseline (Hormozdiari et al. 2018). We calculated coefficient P-values from the coefficient z-scores using a one-sided test assuming a standard normal distribution. We calculated FDR separately for enrichment P-values and coefficient P-values using the BH procedure (Benjamini and Hochberg 1995) and considered traits with FDR < 5% as significantly enriched.

Colocalization of GWAS signals with caQTLs and eQTLs

We downloaded GWAS summary statistics and signal lead variants for 17 traits: three liver enzyme traits (Pazoki et al. 2021), four lipid/cholesterol traits (Graham et al. 2021), four glycemic traits (Chen et al. 2021), type 2 diabetes (Mahajan et al. 2018), BMI (Yengo et al. 2018), waist–hip ratio adjusted for BMI (Pulit et al. 2019), vitamin D (Revez et al. 2020), coronary artery disease (Aragam et al. 2022), and CRP (Said et al. 2022). Although some of these GWAS studies consisted of multiple populations, we used lead variants and summary statistics from individuals genetically similar to European populations for all GWAS studies. We used conditionally distinct signal lead variants if they were provided by the study. For the two studies that did not report conditionally distinct signals (Graham et al. 2021; Said et al. 2022), we identified conditionally distinct lead variants using the following approach. We first defined loci from the GWAS marginal summary statistics by selecting the most significant variant with P < 5 × 10⁻⁸ in a region and removed all other variants within 500 kb (see Supplemental Methods). Loci whose leads were within 1 Mb of each other were merged. We then identified conditionally distinct lead variants within each locus using GCTA (Yang et al. 2011) cojo-slct using variants with MAF ≥ 1%, a collinearity threshold of 0.5, and LD estimated from about 40,000 UKBB donors genetically similar to European populations and extending 500 kb on either side of the lead variant(s) for the locus. For all studies, we performed colocalizations with conditional summary statistics generated using GCTA (Yang et al. 2011) cojo-cond for all signals in which the lead variant was within 500 kb of the lead variant of another signal; we performed colocalizations using marginal summary statistics otherwise. We used liftOver (Hinrichs et al. 2006) to convert coordinates of GWAS variants to GRCh38 to match the caQTL and eQTL data. For caQTLs, we used the summary statistics from all donors. For eQTLs, we used the summary statistics from donors similar to European populations and restricted to protein-coding and lncRNA genes. We tested for colocalization between a GWAS signal and either a caQTL or eQTL signal using coloc.abf (Giambartolomei et al. 2014) if the lead variants of the signals had LD r² ≥ 0.5 (TOPMed Europeans). We ran coloc using betas and standard errors, and only used variants with MAF ≥ 1%. We considered signals that had posterior probability of colocalization (PP H4) > 0.8 to be colocalized. We used LocusZoom v1.4 (Pruim et al. 2010) to make plots of colocalized signals.

Because GWAS signals may be shared across traits, we estimated the number of unique colocalized GWAS signals across traits by first LD-clumping (r² > 0.5, 1000G Europeans) (The 1000 Genomes Project Consortium 2015) GWAS leads across all traits and then counting the number of clumped signals that had at least one variant colocalized with a caQTL. We used the r² > 0.5 threshold to avoid overcounting signals.

Prioritizing caQTLs signals for functional follow-up

We considered caQTLs for functional follow-up if they colocalized with an eQTL and a GWAS signal for lipid or liver enzyme traits and if the caPeak was in a set of coordinately regulated peaks. Because our CRISPRi experiments were performed in HepG2, we additionally required that the caPeak overlap an HepG2 ATAC peak (ENCODE accession ENCSR291GJU) (The ENCODE Project Consortium et al. 2020), that the predicted target gene from the eQTL was expressed in HepG2 (ENCODE accession ENCSR181ZGR), and that the genotype of the caQTL lead variant in HepG2 (ENCODE accession ENCFF989GDN) was either heterozygous or homozygous for the more accessible allele (Zhou et al. 2019).

Transcriptional reporter luciferase assays

We tested for allelic differences in transcriptional activity in HepG2 cells (López-Terrada et al. 2009) for variants overlapping caPeaks at the TENM2 and RALGPS2 loci. HepG2 cells were cultured in MEM-alpha supplemented with 10% FBS and 1 mM sodium pyruvate and maintained at 37°C with 5% CO₂. We designed PCR primers to amplify DNA elements spanning predicted regulatory variants within caPeaks (Supplemental Table 25). For RALGPS2, we generated PCR products for both alleles using DNA from individuals homozygous or heterozygous for the variants, and for TENM2, we also used site-directed mutagenesis (Agilent QuikChange) to generate products with the alternate allele and to make other variable sites identical. We cloned the products into luciferase reporter vector pGL4.23 (Promega) as previously described (Fogarty et al. 2014). The day before transfection, we plated 90,000–120,000 cells into collagen-coated 24-well plates. We transfected 8 wells for each of four sequence-verified constructs for TENM2, and duplicate wells with eight sequence-verified constructs for RALGPS2. We cotransfected wells with phRL-TK Renilla reporter vector using lipofectamine 3000 (Life Technologies). We measured firefly luciferase activity 48 h after transfection and normalized to Renilla activity. We calculated the FC in luciferase activity relative to an empty pGL4.23 vector and used two-tailed t-tests to detect allelic differences in activity. We repeated experiments on a separate day and obtained similar results.

Selection and cloning of gRNA lentiviruses

We designed 10 forward and reverse gRNAs to target a 430 bp region (Chr 1: 178,551,259–178,551,689) within caPeak21014 using CRISPOR (Supplemental Table 26; http://crispor.tefor.net/). We used nontargeting control (NTC) gRNAs as described previously (Smith et al. 2022). We cloned all gRNAs as a pool into the Lenti U6-sgRNA/EF1a-mCherry vector (a gift from Jeremy Day, Addgene 114199) using BbsI-HF and a modified cloning protocol (Ran et al. 2013). Briefly, we annealed 10 µM forward and reverse gRNA primer pairs, phosphorylated using T4 PNK and T4 ligase buffer (NEB), and pooled the 10 gRNA blocks. We ligated 30 ng of BbsSI-HF-digested and calf intestinal phosphatase (CIP)–treated vector with 0.05 µM of annealed and phosphorylated gRNA oligos using T4-DNA ligase (NEB) for 4 h at 16°C, used 2 µL of the ligated product to transform 25 µL of stable competent cells (NEB), plated 125 µL of transformed cells onto ampicillin agar plates, gently scraped all more than 1000 colonies into 250 mL LB broth with ampicillin, and performed overnignt shaking at 37°C. We prepared plasmid DNA (HiPure plasmid maxiprep kit, Invitrogen) and confirmed the presence of all gRNAs by colony PCR.

We produced gRNA lentiviruses as described previously (Perrin et al. 2021). Briefly, we grew 4 × 10⁶ HEK293T cells on a 10 cm plate; cotransfected with 9.5 µg lenti U6-sgRNA/EF1a-mCherry construct, 8 µg of packaging plasmid (psPAX2, a gift from Didier Trono, Addgene plasmid 12260), and 2.5 µg of an envelope plasmid (pMD2.G, a gift from Didier Trono, Addgene plasmid 12259) using Lipofectamine 2000 and PLUS reagent (Invitrogen); and replaced growth media after 18 h. We harvested viral supernatants at 48 and 72 h after transfection and concentrated using Lenti-X concentrators (Clontech). We functionally tittered the lentivirus by amplifying viral DNA WPRE, plasmid backbone DNA, and genomic DNA (LP34) for each sample as previously described (Gordon et al. 2020) and represented functional titers as multiplicity of infection (MOI).

CRISPRi

We used a doxycycline-inducible CRISPRi system (López-Terrada et al. 2009) by stably expressing dCas9 fused to a KRAB-repressor domain in the AAVS1 safe-harbor locus (Hazelbaker et al. 2020) as previously described (Pandey et al. 2024). We induced HepG2-dCas9-KRAB cells with doxycycline (2 µg/mL) for 48 h, sorted for GFP using a Becton Dickinson FACSAria II at the UNC Flow Cytometry Core Facility, and plated 1.2 × 10⁵ GFP-positive cells/well onto collagen-I (Corning)-coated 12-well plates. We maintained cells in 275 µg/mL neomycin except during the transduction period. We spin-infected cells with lenti-gRNAs (Perrin et al. 2021) at 40 MOI in the presence of 10 µg/mL polybrene. After 8 h, we replaced growth media with media containing doxycycline (2 µg/mL) and neomycin and visualized GFP and mCherry expression by fluorescent microscopy. Ninety-six hours after transduction, we lysed cells, extracted total RNA (RNeasy plus, Qiagen), and converted to cDNA (Superscript IV Vilo, Invitrogen). We assessed gene expression by quantitative PCR using TaqMan probes for RALGPS2 (Thermo Fisher Scientific) and normalized to B2M (Hs00187842-m1) expression level. We performed experiments on two separate days and obtained similar results. For the second experiment, we removed one replicate from the nontargeting control experiments that was more than two standard deviations from the mean of gene expression.