Published papers

Objectives Synovial fibroblasts (SFs) are one of the major components of the inflamed synovium in rheumatoid arthritis (RA). We aimed to gain insight into the pathogenic mechanisms of SFs through elucidating the genetic contribution to molecular regulatory networks under inflammatory condition. Methods SFs from RA and osteoarthritis (OA) patients (n=30 each) were stimulated with eight different cytokines (interferon (IFN)-α, IFN- γ, tumour necrosis factor-α, interleukin (IL)-1β, IL-6/sIL-6R, IL-17, transforming growth factor-β1, IL-18) or a combination of all 8 (8-mix). Peripheral blood mononuclear cells were fractioned into five immune cell subsets (CD4 + T cells, CD8 + T cells, B cells, natural killer (NK) cells, monocytes). Integrative analyses including mRNA expression, histone modifications (H3K27ac, H3K4me1, H3K4me3), three-dimensional (3D) genome architecture and genetic variations of single nucleotide polymorphisms (SNPs) were performed. Results Unstimulated RASFs differed markedly from OASFs in the transcriptome and epigenome. Meanwhile, most of the responses to stimulations were shared between the diseases. Activated SFs expressed pathogenic genes, including CD40 whose induction by IFN-γwas significantly affected by an RA risk SNP (rs6074022). On chromatin remodelling in activated SFs, RA risk loci were enriched in clusters of enhancers (super-enhancers; SEs) induced by synergistic proinflammatory cytokines. An RA risk SNP (rs28411362), located in an SE under synergistically acting cytokines, formed 3D contact with the promoter of metal-regulatory transcription factor-1 (MTF1) gene, whose binding motif showed significant enrichment in stimulation specific-SEs. Consistently, inhibition of MTF1 suppressed cytokine and chemokine production from SFs and ameliorated mice model of arthritis. Conclusions Our findings established the dynamic landscape of activated SFs and yielded potential therapeutic targets associated with genetic risk of RA.

Objective: Considering the diverse aetiologies and immunodysregulatory statuses observed in each patient with rheumatoid arthritis (RA), stratification based on peripheral blood immunophenotyping holds the potential to enhance therapeutic responses to molecular targeted therapies, biological/targeted synthetic disease-modifying antirheumatic drugs (b/tsDMARDs). Methods: Immunophenotype analysis was conducted on a cohort of over 500 b/tsDMARDs-naïve patients using flow cytometry. Patients with RA were stratified based on their immunophenotypes, and the treatment response to each targeted therapy was evaluated. Validation was performed using an additional cohort of 183 b/tsDMARDs-naïve patients with RA. Results: Patients with RA were stratified into five clusters, two of which exhibited distinct RA phenotypes compared with controls, characterised by significant increases in CD4+ effector memory T cells re-expressing CD45RA. Notably, the effectiveness of different b/tsDMARDs varied across clusters. The group using promising b/tsDMARDs was labelled as ‘expected’ whereas the ‘non-expected’ group comprised those using others. The expected group outperformed the non-expected group with higher 26-week remission rates (39.9% vs 24.6%, p=0.0004) and low disease activity achievement (80.8% vs 60.2%, p<0.0001). Trajectory analysis showed the non-expected group's 26-week disease activity was influenced by Clinical Disease Activity Index at baseline unlike the expected group. Additionally, different molecular targeted therapies influenced the proportions of each immune cell subset variably. To validate, immunophenotyping was performed on a validation cohort. When 183 cases were grouped based on their b/tsDMARDs usage into expected/non-expected groups, the expected group had a higher remission rate (p=0.0021), further confirming the observed trend. Conclusion: Our findings offer valuable insights into the diversity of RA and potential therapeutic strategies grounded in the molecular underpinnings.

To improve imputation quality for genome-wide association studies (GWAS) conducted on the Japanese population, we developed and evaluated four Japanese population-specific reference panels. These panels were constructed through the augmentation of the 1000 Genomes Project (1KG) panel using Japanese whole genome sequencing (WGS) data, with sample sizes ranging from 1 K to 7 K individuals enrolled through the Biobank Japan (BBJ) project, and sequencing depths ranging from 3× to 30×. Among these panels, an augmented reference panel comprising 7472 WGS samples of mixed depth (1KG+7K) exhibit the greatest improvement in imputation quality relative to the Trans-Omics for Precision Medicine (TOPMed) reference panel. Notably, we observe these improvements primarily for rare variants with a minor allele frequency (MAF) <5%. To demonstrate the benefits of improved imputation quality in association analyses of complex traits, we conducted GWAS for serum uric acid and total cholesterol levels following imputation up to the 1KG+7K panel. The analysis reveals several loci reaching genome-wide significance (P < 5 × 10–8) in the 1KG+7K imputation output yet remaining undetected when the same sample set is imputed up to the TOPMed reference panel. In summary, the 1KG+7K panel demonstrates significant advantages in the discovery of trait-associated loci, particularly those influenced by low-frequency association signals.

Gene transcription is regulated through complex mechanisms involving non-coding RNAs (ncRNAs). As the transcription of ncRNAs, especially of enhancer RNAs, is often low and cell type specific, how the levels of RNA transcription depend on genotype remains largely unexplored. Here we report the development and utility of a machine-learning model (MENTR) that reliably links genome sequence and ncRNA expression at the cell type level. Effects on ncRNA transcription predicted by the model were concordant with estimates from published studies in a cell-type-dependent manner, regardless of allele frequency and genetic linkage. Among 41,223 variants from genome-wide association studies, the model identified 7,775 enhancer RNAs and 3,548 long ncRNAs causally associated with complex traits across 348 major human primary cells and tissues, such as rare variants plausibly altering the transcription of enhancer RNAs to influence the risks of Crohn’s disease and asthma. The model may aid the discovery of causal variants and the generation of testable hypotheses for biological mechanisms driving complex traits.

T cells acquire a regulatory phenotype when their T cell antigen receptors (TCRs) experience an intermediate- to high-affinity interaction with a self-peptide presented via the major histocompatibility complex (MHC). Using TCRβ sequences from flow-sorted human cells, we identified TCR features that promote regulatory T cell (Treg) fate. From these results, we developed a scoring system to quantify TCR-intrinsic regulatory potential (TiRP). When applied to the tumor microenvironment, TiRP scoring helped to explain why only some T cell clones maintained the conventional T cell (Tconv) phenotype through expansion. To elucidate drivers of these predictive TCR features, we then examined the two elements of the Treg TCR ligand separately: the self-peptide and the human MHC class II molecule. These analyses revealed that hydrophobicity in the third complementarity-determining region (CDR3β) of the TCR promotes reactivity to self-peptides, while TCR variable gene (TRBV gene) usage shapes the TCR’s general propensity for human MHC class II-restricted activation.

Non-coding genetic variants may cause disease by modulating gene expression. However, identifying these expression quantitative trait loci (eQTLs) is complicated by differences in gene regulation across fluid functional cell states within cell types. These states—for example, neurotransmitter-driven programs in astrocytes or perivascular fibroblast differentiation—are obscured in eQTL studies that aggregate cells1,2. Here we modelled eQTLs at single-cell resolution in one complex cell type: memory T cells. Using more than 500,000 unstimulated memory T cells from 259 Peruvian individuals, we show that around one-third of 6,511 cis-eQTLs had effects that were mediated by continuous multimodally defined cell states, such as cytotoxicity and regulatory capacity. In some loci, independent eQTL variants had opposing cell-state relationships. Autoimmune variants were enriched in cell-state-dependent eQTLs, including risk variants for rheumatoid arthritis near ORMDL3 and CTLA4; this indicates that cell-state context is crucial to understanding potential eQTL pathogenicity. Moreover, continuous cell states explained more variation in eQTLs than did conventional discrete categories, such as CD4+ versus CD8+, suggesting that modelling eQTLs and cell states at single-cell resolution can expand insight into gene regulation in functionally heterogeneous cell types.

Alternative splicing contributes to complex traits, but whether this differs in trait-relevant cell types across diverse genetic ancestries is unclear. Here we describe cell-type-specific, sex-biased and ancestry-biased alternative splicing in ~1 M peripheral blood mononuclear cells from 474 healthy donors from the Asian Immune Diversity Atlas. We identify widespread sex-biased and ancestry-biased differential splicing, most of which is cell-type-specific. We identify 11,577 independent cis-splicing quantitative trait loci (sQTLs), 607 trans-sGenes and 107 dynamic sQTLs. Colocalization between cis-eQTLs and trans-sQTLs revealed a cell-type-specific regulatory relationship between HNRNPLL and PTPRC. We observed an enrichment of cis-sQTL effects in autoimmune and inflammatory disease heritability. Specifically, we functionally validated an Asian-specific sQTL disrupting the 5′ splice site of TCHP exon 4 that putatively modulates the risk of Graves’ disease in East Asian populations. Our work highlights the impact of ancestral diversity on splicing and provides a roadmap to dissect its role in complex diseases at single-cell resolution.

Splicing quantitative trait loci (sQTLs) are one of the major causal mechanisms in genome-wide association study (GWAS) loci, but their role in disease pathogenesis is poorly understood. One reason is the complexity of alternative splicing events producing many unknown isoforms. Here, we propose two approaches, namely integration and selection, for this complexity by focusing on protein-structure of isoforms. First, we integrate isoforms with the same coding sequence (CDS) and identify 369-601 integrated-isoform ratio QTLs (i2-rQTLs), which altered protein-structure, in six immune subsets. Second, we select CDS incomplete isoforms annotated in GENCODE and identify 175-337 isoform-ratio QTL (i-rQTL). By comprehensive long-read capture RNA-sequencing among these incomplete isoforms, we reveal 29 full-length isoforms with unannotated CDSs associated with GWAS traits. Furthermore, we show that disease-causal sQTL genes can be identified by evaluating their trans-eQTL effects. Our approaches highlight the understudied role of protein-altering sQTLs and are broadly applicable to other tissues and diseases.

Synovial tissue inflammation is a hallmark of rheumatoid arthritis (RA). Recent work has identified prominent pathogenic cell states in inflamed RA synovial tissue, such as T peripheral helper cells; however, the epigenetic regulation of these states has yet to be defined. Here, we examine genome-wide open chromatin at single-cell resolution in 30 synovial tissue samples, including 12 samples with transcriptional data in multimodal experiments. We identify 24 chromatin classes and predict their associated transcription factors, including a CD8 + GZMK+ class associated with EOMES and a lining fibroblast class associated with AP-1. By integrating with an RA tissue transcriptional atlas, we propose that these chromatin classes represent ‘superstates’ corresponding to multiple transcriptional cell states. Finally, we demonstrate the utility of this RA tissue chromatin atlas through the associations between disease phenotypes and chromatin class abundance, as well as the nomination of classes mediating the effects of putatively causal RA genetic variants.

Autoimmune disease heritability is enriched in T cell-specific regulatory regions of the genome. Modern-day T cell datasets now enable association studies between single nucleotide polymorphisms (SNPs) and a myriad of molecular phenotypes, including chromatin accessibility, gene expression, transcriptional programs, T cell antigen receptor (TCR) amino acid usage, and cell state abundances. Such studies have identified hundreds of quantitative trait loci (QTLs) in T cells that colocalize with genetic risk for autoimmune disease. The key challenge facing immunologists today lies in synthesizing these results toward a unified understanding of the autoimmune T cell: which genes, cell states, and antigens drive tissue destruction?

Translating genome-wide association study (GWAS) loci into causal variants and genes requires accurate cell-type-specific enhancer–gene maps from disease-relevant tissues. Building enhancer–gene maps is essential but challenging with current experimental methods in primary human tissues. Here we developed a nonparametric statistical method, SCENT (single-cell enhancer target gene mapping), that models association between enhancer chromatin accessibility and gene expression in single-cell or nucleus multimodal RNA sequencing and ATAC sequencing data. We applied SCENT to 9 multimodal datasets including >120,000 single cells or nuclei and created 23 cell-type-specific enhancer–gene maps. These maps were highly enriched for causal variants in expression quantitative loci and GWAS for 1,143 diseases and traits. We identified likely causal genes for both common and rare diseases and linked somatic mutation hotspots to target genes. We demonstrate that application of SCENT to multimodal data from disease-relevant human tissue enables the scalable construction of accurate cell-type-specific enhancer–gene maps, essential for defining noncoding variant function.

The fluctuating nature of disease activity in systemic lupus erythematosus (SLE), alternating between flares and remissions, poses substantial challenges for its effective management. The use of current biomarkers for monitoring SLE is limited in clinical settings owing to insufficient comprehension of the complex immune involvement underlying the disease course. Here, therefore, we profiled peripheral blood mononuclear cells at both stable and exacerbation states (total of n = 19) from six patients with SLE and 32 healthy donors using integrated single-cell RNA and T cell receptor (TCR) sequencing. To validate our findings, we analyzed two independent external datasets: bulk RNA sequencing and TCR data from 79 controls and 62 patients with SLE and single-cell RNA sequencing data from 99 healthy controls and 162 patients with SLE. Our analysis revealed cell type-specific activation of interferon-related genes in SLE grouped into four clusters, with elevated activity in disease-associated immune cells. Among these, atypical B cells associated with autoantibody production exhibited distinct differentiation patterns compared with conventional memory B cells, driven by heightened interferon signaling in SLE. Notably, clonal expansion of effector CD8 T cells emerged as a key driver of disease exacerbation, as indicated by reduced TCR diversity. Specific CD8 T cell clonotypes expanded during flare states, transitioning to effector phenotypes that exhibited heightened cytotoxicity and amplified interferon signaling, strongly correlating with tissue damage and flare severity. Our findings establish a critical link between interferon-driven mechanisms and cytotoxic T cell dysfunction in SLE flares, offering potential targets for therapeutic intervention and predictive biomarkers.