Single-nucleus chromatin and gene expression profiling across hundreds of skeletal muscle samples reveals context-specific regulation

Skeletal muscle is relevant to several polygenic metabolic traits and diseases including type 2 diabetes (T2D). Identifying genetic mechanisms underlying these traits requires pinpointing the relevant cell-types, regulatory elements, target genes, and causal variants. Here, we used genetic multiplexing to generate population-scale single nucleus (sn) chromatin accessibility (snATAC-seq) and transcriptome (snRNA-seq) maps across 287 frozen human skeletal muscle biopsies representing 456,880 nuclei. We identified 13 cell-types that collectively represented 983,155 ATAC summits. We integrated genetic information to identify 6,866 expression quantitative trait loci (eQTL) and 100,928 chromatin accessibility QTL (caQTL) (5% FDR) across the five most abundant cell-types, and show atlas-level snATAC maps often miss peaks. We identified 12,046 e-caQTL colocalizations and inferred causality between chromatin accessibility and gene expression through mediation analyses. 3,378 genome-wide association study (GWAS) signals across 43 relevant traits colocalized with sn e/caQTL, 52% in a cell-type-specific manner. 77% (3,616) GWAS signals colocalized with caQTLs and not eQTL, highlighting the importance of sn-caQTL maps for GWAS functional studies. GWAS-e/caQTL colocalization showed distinct cell-type-specific regulatory paradigms. For example, a C2CD4A/B T2D GWAS signal colocalized with caQTLs in muscle fibers and imputed chromatin contacts nominated VPS13C - a glucose uptake gene. The caQTL peak overlapping caSNP rs7163757 showed allelic regulatory activity differences in a human myocyte cell line massively parallel reporter assay. These results illuminate the genetic regulatory architecture of human skeletal muscle at high-resolution epigenomic, transcriptomic, and cell state scales and serve as a template for population-scale multi-omic mapping in complex tissues and traits. Overall design: Massively Parallel Reporter Assay experimenatal design: Cloning: We ordered oligos as 230 bp sequences where 197 bp comprise the variant of interest flanked on both by 98 bp of genomic context, and the additional 33 bp are cloning adapters. Within this panel, we included a set of ~50 negative control sequences defined by a previous publication107 We added 20 bp barcodes via a 2-step PCR amplification process then incorporated the barcoded oligos into a modified pMPRA1 vector (a gift from Tarjei Mikkelsen108, Addgene #49349) upstream of the GFP reporter gene using Golden Gate assembly. After transforming and expanding in NEB 10-beta electrocompetent bacteria, we sequenced this version of the MPRA library to establish a barcode-oligo pairing dictionary. We performed a second Golden Gate assembly step to insert an ENCODE-annotated promoter for the human MYBPC2 gene in between the oligo and barcode. Finally, we used restriction cloning to port the assembled MPRA block (oligo, barcode, promoter, GFP) to a lentiviral transfer vector, which was used by the University of Michigan viral vector core to produce infectious lentiviral particles. Primer sequences used for cloning and sequencing library preparation along with the MYBPC2 promoter sequence are included in a separate table. MPRA Experiment For each replicate, we infected 4x106 LHCN-M2 human skeletal myoblasts with our MPRA library at an MOI of ~10. After infection, we passaged the cells for one week to remove any unincorporated virus or contaminating transfer plasmid, then differentiated the cells for one week. We isolated RNA and gDNA from each replicate using the Qiagen AllPrep DNA/RNA mini kit. We reverse transcribed RNA into cDNA with a GFP-specific primer, then constructed indexed sequencing libraries for both the cDNA and gDNA libraries using Illumina-compatible primers.859 Data Analysis After quality checks and filtering, we calculated the sum of barcode counts for each oligo within a replicate. We used DESeq2 v1.34.094 to perform normalization and differential expression analysis. We used a nested model to identify oligos with significant activity (relative to plasmid input) and significant allelic bias (between reference and alternate alleles). All results were subject to a Benjamini-Hochberg FDR of 5%