that the eQTLs result from independent var-
iants that influence expression in different cell
types. To test between these hypotheses, we
performed a regression strategy to evaluate the
change in the test statistic of an eSNP after
regressing out the effects of an eSNP from
another cell type. Under this strategy, if the
eSNPs tag the same causal variant for that
gene or are in linkage disequilibrium with one
another, then the allelic effect size of the orig-
inal eSNP will decrease in the conditional analy-
sis. Similarly, if they tag independent variants,
the allelic effect will remain relatively unchanged.
We performed this strategy for each pairwise
combination of eQTLs where different top eSNPs
were identified in different cell types.
We tested whether each eGene was tagged
by two distinct variants by conditioning the
lead eSNP from the first cell type on the lead
eSNP from the second cell type for every pair
of cell types (182 pairs). The correlation coef-
ficients of significant independent eSNPs from
shared eGenes pre- and postconditioning are
shown in Fig. 2D and fig. S16. Whereas most
lymphoid immune cell eQTLs had a consider-
able change in the correlation coefficients
after conditioning, among the myeloid im-
mune cells, the eQTL correlation coefficients
remained similar (fig. S16). This finding sug-
gests that lymphoid cell types are more likely
to share genetic control of gene expression be-
tween cell types compared with myeloid cells.
Evidence suggests that cell type–specific
chromatin accessibility underlies a proportion
of cell type–specific cis-eQTLs
To explore the functional regulation underly-
ing cis-eQTLs, we tested for the overlap of
eSNP locations and regions of open chromatin
generated from single-cell assay for transpos-
ase accessible chromatin sequencing (scATAC-
seq) data from 8876 cells. Cells were classified
into each of the 14 cell types, and open chro-
matin peaks were called for each cell type that
had more than five classified cells. This filter-
ing retained 11 cell types, comprising the most
abundant populations [except CD4+T cells with
an effector memory or central memory pheno-
type (CD4ET), CD4SOX4, and plasma cells] (fig.
S18). On average, we identified 52,048 peaks
per cell type, with the mean distance between
an eSNP and the nearest peak ranging from
7485 to 31,383 base pairs. To determine whether
the location of cis-eQTLs was significantly closer
to open chromatin regions, we compared the
distances between the cis-eQTLs. We random-
ly sampled SNPs that were selected based on
the same distance distribution from the tran-
script to the nearest peaks per cell type using
a bootstrapping technique. We observed a
significant difference between the cis-eQTL
distances across all cell types except CD4SOX4
cells [false discovery rate (FDR) < 0.05] (fig.
S19). We conclude from these results that cell
type–specific chromatin accessibility is likely
to contribute to variation in allelic effects on
gene expression between cell types.Single-cell eQTLs replicate in multiethnic
cohorts and bulk eQTL studies
To verify cell-specific eQTL findings, we repli-
cated our lead eSNP results in two indepen-
dent cohorts of European and Asian ancestry,
consisting of 113 and 89 individuals, respec-
tively. Of the 16,597 eSNP 1 -eGene pairs, 10,071
were present with a minor allele frequency
greater than 0.05 in both cohorts. Of these,
3198 (26%) in the European cohort and 2243
(22%)intheAsiancohortreplicatedattheFDR
threshold of 5%, which is encouraging given the
differences between the sample sizes of these
cohorts and the sample size of the OneK1K
discovery cohort (tables S5, S12, and S13).
Indeed, correcting the FDR distributions
under the assumptions of equal sample size in
the discovery and replication cohorts leads to
87 and 78% replication rates in the European
and Asian cohorts, respectively. Similarly, the
concordance of allelic direction over all tested
loci was 76.0 to 98.1% in the European cohort
and 72.2 to 95.4% in the Asian cohort. This
concordance increases to 99.3 to 100% and
96.9 to 99.8%, respectively, for eQTLs replicat-
ing at an FDR less than 0.05 (Fig. 3, A and B).
The discrepancy in replication rates between
cohorts likely reflects differences in the allele
frequencies of eSNPs between population
groups. However, the results indicate that cell
type–specific eQTLs are likely to be largely
shared among populations. The discovery of
OneK1K eQTLs was tested for replication in all
cell types in the replication cohorts. At an FDR
less than 0.05, replicating eQTLs and eGenes
are predominantly identified in a single cell
type (Fig. 3, C and D), providing further evi-
dence for cell type–specific effects of loci on
gene expression in PBMCs. The concordance of
correlation coefficients between the OneK1K
and replication cohorts are shown in Fig. 3E for
both the European and Asian samples. We
were able to replicate 62.5 and 40.4% of cis-
eQTLs identified in bulk RNA-seq studies of
blood samples from the eQTL-Gen Consortium
( 34 ) and GTEx Consortium ( 5 ), respectively
(fig. S20 and table S14).Identification of dynamic eQTL allelic effects
across the B cell landscape
We investigated the dynamic effects of eQTLs
across the pseudotime landscape of immature
and naïve B (BIN) cells through to memory B
(BMem) cells. Cells were categorized into six
quantiles (Q1 to Q6) based on their relative
position on the pseudotime curve (Fig. 4, A
and B). Overlaying the expression of classical
markers revealed a graded change across the
derived trajectory from BIN(Q1) to BMemcells
(Q6). For example,TCL1AandIL4Rare highlyexpressed in naïve B cells ( 35 , 36 ) and were
found to be down-regulated across the transi-
tion to BMemcells (Fig. 4C). Conversely, the
expression ofCD27, a canonical BMemcell
marker ( 37 ), increased as the cells transitioned
to a memory state.IgJexpression, a marker of
immunoglobulin M (IgM) and IgA production,
was up-regulated in the higher quantiles, in-
dicating that they contain cells poised to be-
come plasma cells ( 38 ) (Fig. 4C).
We sought to identify instances where eQTL
allelic effects exhibited either linear or nonlinear
changes across the trajectory of naïve to mem-
ory B cell transition. Dynamic B cell eQTLs were
determined by testing the interaction between
the genotype and quantile ranks using both
linear and quadratic models. Of the 3074 cis-
eQTLs identified in BINand BMemcells, 1988
were expressed in at least three pseudotime
quantiles and tested for dynamic effects. Of
these, we identified significant changes in the
allelic effect across the trajectory for 333 of
them (FDR < 0.05) (Fig. 4D and fig. S21).
Many of the genes with dynamic eQTL ef-
fects have a role in fine-tuning B cell migration,
activation, survival, or function. For example,
SELLis involved in integrin-mediated migra-
tion to and within tissues ( 39 , 40 ). Migration
to and organization of B cells within the germ-
inal center is a critical component in generating
appropriate memory and humoral outputs. The
allelic effect of the intronic variant rs4987360-G
onSELLexpression is largest in immature
cells, decreasing over each of the subsequent
quantiles (Fig. 4E). The opposite trend is iden-
tified for SNPs that influence the expression
of the Src family tyrosine kinase B lymphocyte
kinase (BLK), a gene responsible for regulating
the amplitude of signaling downstream of the
B cell receptor. Both rs2736336 and rs2409780
show the greatest allelic effects in Q5 and Q6
(Fig. 4E and table S15). Interestingly, rs2736336,
a variant in the promoter ofBLK, is associated
with systemic lupus erythematosus (SLE) ( 41 ),
whereas rs2409780, an intronic variant, is in
high linkage disequilibrium with variants
associated with SLE and RA [coefficient of
determination (R^2 ) = 0.99, and coefficient of
linkage disequilibrium (D′)=0.99]( 13 , 42 ).
Another gene responsible for interpreting
the signaling downstream of B cell surface
receptors and influencing subsequent B cell
proliferation and survival is c-Rel, encoded by
the transcription factorREL( 43 ). rs12989427
is in high linkage disequilibrium with variants
associated with SLE (R^2 = 0.88, andD′= 0.98),
and the allelic effect follows a nonlinear rela-
tionship, peaking at the medium point of the B
cell trajectory (Fig. 4E).ORMDL3promotes
mature B cell survival by suppressing apopto-
sis and promoting autophagy ( 44 ). rs7359623
and rs8067378 are in high linkage disequil-
ibrium with risk variants (R^2 > 0.8, andD′>
0.9) implicated in a range of autoimmuneYazaret al.,Science 376 , eabf3041 (2022) 8 April 2022 5 of 14
RESEARCH | RESEARCH ARTICLE