Article Text
PDF
Abstract
Background Although genome-wide association studies (GWASs) have identified more than 100 loci associated with rheumatoid arthritis (RA) susceptibility, the causal genes and biological mechanisms remain largely unknown.
Methods A cross-tissue transcriptome-wide association study (TWAS) using the unified test for molecular signaturestool was performed to integrate GWAS summary statistics from 58 284 individuals (14 361 RA cases and 43 923 controls) with gene-expression matrix in the Genotype-Tissue Expression project. Subsequently, a single tissue by using FUSION software was conducted to validate the significant associations. We also compared the TWAS with different gene-based methodologies, including Summary Data Based Mendelian Randomization (SMR) and Multimarker Analysis of Genomic Annotation (MAGMA). Further in silico analyses (conditional and joint analysis, differential expression analysis and gene-set enrichment analysis) were used to deepen our understanding of genetic architecture and comorbidity aetiology of RA.
Results We identified a total of 47 significant candidate genes for RA in both cross-tissue and single-tissue test after multiple testing correction, of which 40 TWAS-identified genes were verified by SMR or MAGMA. Among them, 13 genes were situated outside of previously reported significant loci by RA GWAS. Both TWAS-based and MAGMA-based enrichment analyses illustrated the shared genetic determinants among autoimmune thyroid disease, asthma, type I diabetes mellitus and RA.
Conclusion Our study unveils 13 new candidate genes whose predicted expression is associated with risk of RA, providing new insights into the underlying genetic architecture of RA.
- Rheumatoid Arthritis
- Polymorphism, Genetic
- Epidemiology
Data availability statement
The datasets used and/or analysed during the current study are available from the corresponding authors on reasonable request.
This is an open access article distributed in accordance with the Creative Commons Attribution 4.0 Unported (CC BY 4.0) license, which permits others to copy, redistribute, remix, transform and build upon this work for any purpose, provided the original work is properly cited, a link to the licence is given, and indication of whether changes were made. See: https://creativecommons.org/licenses/by/4.0/.
Statistics from Altmetric.com
WHAT IS ALREADY KNOWN ON THIS TOPIC
Although genome-wide association studies (GWASs) have identified more than 100 rheumatoid arthritis (RA) risk loci, these variants can only explain 15% of RA heritability.
Most GWAS signals are located in non-coding regions, where often overlap with gene regulatory elements and highly enrich expression quantitative trait loci.
WHAT THIS STUDY ADDS
Through a two-stage transcriptome-wide association study design, we unveil 13 novel susceptibility genes whose genetically predicted expression is associated with risk of RA.
Our study reveals shared genetic determinants among autoimmune thyroid disease, asthma, type I diabetes mellitus and RA.
HOW THIS STUDY MIGHT AFFECT RESEARCH, PRACTICE OR POLICY
Our finding about novel RA loci refines the known genetic architecture of RA and provides genetic biomarkers for RA.
Introduction
Rheumatoid arthritis (RA), a chronic autoimmune disease, is characterised by persistent synovial inflammation affecting around 0.5%–1.0% of general population worldwide.1 Currently, there is no completely curative therapy available for this highly disabling disease. Hence, a better knowledge of the underlying mechanisms of RA will contribute to the development of effective therapeutic targets.
Based on twin studies, the relative contribution of genetic variation to the liability of having RA has been estimated to be around 60%.2 Although genome-wide association studies (GWASs) have identified more than 100 RA risk loci, these variants can only explain 15% of RA heritability.3 In addition, the biological interpretation and functional comprehension of these associations remain elusive.4 5 Most GWAS signals are located in non-coding regions, where often overlap with gene regulatory elements and highly enrich expression quantitative trait loci (eQTL).6 7 These clues point to a crucial role of transcriptional regulation in impacting RA susceptibility.
Transcriptome-wide association study (TWAS) leverages eQTL and GWAS data to identify novel susceptibility genes whose genetically predicted expression is associated with disease risk.8 This method aggregates several variants into a functional gene unit, reducing the number of multiple comparisons and mapping to the candidate gene directly. As a powerful gene-based method, TWAS has been successfully applied to decipher the genetic architecture of several complex traits.9–12 The majority of existing TWAS studies calculated genetic–expression matrix in each tissue separately, which might neglect the sharing local regulation of gene expression across different tissues.13 Evidence has shown that eQTL with large effects can regulate gene expression across multiple tissues.14 A recent cross-tissue TWAS approach, called the unified test for molecular signatures (UTMOST), was developed to perform gene-level association analysis across different tissues, improving the accuracy and effectiveness of imputation model than that of single tissue method.15 By imposing a group-lasso penalty on the effect sizes of single nucleotide polymorphisms (SNPs) across tissues, UTMOST encourages eQTL that are shared across tissues and keeps tissue-specific eQTL with strong effects. In recent years, cross-tissue association analysis have been widely used in screening candidate susceptibility genes for complex multisystemic disorders whose biologically relevant tissues are not entirely clear, such as inflammatory bowel diseases, schizophrenia and Alzheimer’s disease.16 17
In this work, we carried out a cross-tissue TWAS of RA by integrating eQTL data from Genotype-Tissue Expression (GTEx) project with largest European RA GWAS. A single tissue TWAS using another Functional Summary-based Imputation (FUSION) method was adopted to validate candidate susceptibility genes synchronously. Follow-up bioinformatic analyses were performed to explore the biological characterisation for these candidate genes.
Methods
RA GWAS data source
We obtained the summary statistics from the largest meta-analysis of RA GWAS, comprising 14 361 cases and 43 923 controls of European ancestry.18 The participants enrolled in the GWAS meta-analysis were obtained from 18 separately studies. All patients were diagnosed with RA by professional rheumatologists according to the 1987 American College of Rheumatology criteria. Detailed description of the quality control, genotyping and imputation procedures were provided in the published studies.19 20
TWAS analyses in cross-tissue and single tissue
A two-stage TWAS design of the study is shown in figure 1. In the discovery stage, we performed a cross-tissue association test by using UTMOST method.15 RA GWAS summary data were integrated with eQTL data of 44 tissues from GTEx to impute the genetic component of gene expression in each tissue. Then, a generalised Berk-Jones test was applied to combine gene–trait associations in 44 GTEx human tissues based on the covariance from single tissue statistics (online supplemental table S1). For each gene, UTMOST trains a cross-tissue expression imputation model based on a penalised multivariate regression, which has considered the different directions and effect sizes of eQTL signals across tissues. To reduce noise in the cross-tissue association test and the risk of false positive rate, we performed a validation to incorporate RA GWAS and eQTL data of whole blood from GTEx by using the FUSION software,8 a widely used method in prior TWAS analysis. FUSION builds predictive models with several penalised linear models (GBLUP, LASSO, Elastic Net, etc) for those significant cis-heritability genes estimated by using SNPs within 500 kb on both sides of the gene boundary, and then chooses the best model based on the R2 calculated by a fivefold cross-validation of each model. In the external validation, we repeated the analysis using the eQTL data of peripheral blood from The Netherlands Twin Register (NTR)21 and calculated Pearson’s correlation coefficient to analyse the effect sizes of testable TWAS genes between GTEx and NTR. TWAS significance for both cross-tissue and single-tissue analysis was established as Benjamini-Hochberg corrected false discovery rate (FDR) value below 0.05.
Supplemental material
Overview of the transcriptome-wide association study design of rheumatoid arthritis (RA). FDR, false discovery rate; GWAS, genome-wide association study; TWAS, transcriptome-wide association study.
Conditional and joint analysis
Conditional and joint analysis for genome-wide FDR-corrected significant TWAS signals was used to evaluate the GWAS association signal after removing the TWAS association signal. Each RA GWAS SNP association was conditioned on the joint gene model. To investigate inflation of imputed association statistics under the null of no GWAS association, the permutation test was employed with a maximum of 100 000 permutations and permutation-based p value threshold of 0.05 for each of the significant loci.
Gene-based association study
Summary data-based Mendelian randomization (SMR) is a novel gene-based approach to prioritise functional genes in GWAS loci, which uses cis-eQTL SNPs, gene expression and RA GWAS as instrumental variables, exposure and outcome, respectively.22 Multimarker analysis of genomic annotation (MAGMA) uses a multiple regression model to calculate the cumulative effect of several SNPs that assigned to a specific gene (±10 kb).23 The phase 3 of 1000 Genomes European population was used for reference panel to calculate the linkage disequilibrium (LD).24
Validation of TWAS results by mRNA expression profiles of RA
To explore whether the expression levels of TWAS identified genes associated with RA risk were dysregulated in RA patients compared with controls, we obtained whole blood gene expression microarray data of RA from Gene Expression Omnibus dataset GSE93272.25 The genome-wide transcript profiling of whole blood was collected from 232 RA patients and 43 healthy controls. The gene expression data were processed with log transformation followed by quantile normalisation using the linear models for microarray data package to identify differentially expressed genes. We selected the differentially expressed genes based on FDR corrected p<0.05.
Pathway enrichment analysis
To uncover the biologically relevant pathways in aetiology of RA, we conducted two gene-set methodologies including TWAS-based gene set enrichment analysis (TWAS-GSEA) and MAGMA-based gene-set analysis synchronously. TWAS-GSEA performs gene enrichment based on a mixed model test, in that it can directly use the results of TWAS. In MAGMA, the second gene-set analysis is established in a linear regression model by using the gene p value and gene correlation matrix. KEGG, BioCarta and Reactome pathways were downloaded from MSigDB database V.7.5.1 (https://www.gsea-msigdb.org/gsea/msigdb).
Results
Transcriptome-wide association study results of RA
For the cross-tissue discovery, a total of 262 genes showed a statistically significant signal after FDR correction (p<0.05) and 96 of them remained significant after Bonferroni correction (p<2.89×10−6) (online supplemental table S2) (online supplemental figure S1). For single-tissue internal validation, of all the 8125 genes modelled in our genotype data with significant cis-heritable expression in the whole blood from GTEx dataset, 187 genes displayed a significant TWAS association signal with PFDR<0.05 (online supplemental table S3). In all, we identified 47 overlapped candidate genes with the stringent threshold in both cross-tissue and single-tissue test, including 13 genes located in novel loci (table 1).
Supplemental material
The significant genes for RA risk in both cross-tissue and single-tissue TWAS analysis
Conditional and joint analysis
As shown in table 2, 11 loci represented the independent signal containing multiple significant genes, including 1p36.32 (TTC34), 5q21.1 (PPIP5K2), 6p21.31 (IP6K3, C6orf106), 6p22.1 (PGBD1, TRIM10), 6p22.2 (HIST1H2BG, ZNF322), 9q33.2 (PHF19, TRAF1), 11q12.2 (TMEM258), 12q13.2 (SUOX), 12q14.1 (TSPAN31), 16p11.2 (SPNS1, RNF40) and 17q21.1 (ORMDL3) (conditional p<0.05). We observed that several GWAS signals were driven by genetically modulated gene expression. For instance, SUOX explained most of the signal at 12q13.2 locus, while TWAS signal for RPS26 dropped obviously when conditioned on the predicted expression of SUOX (figure 2A). Similarly, a cluster of three TWAS significant associations (ORMDL3, GSDMB and PGAP3) were found at 17q21.1 locus, conditional analysis indicated that ORMDL3 explained most of the signals in this region (figure 2B).
Conditional analysis result of multigene loci
Regional association of transcriptome-wide association study (TWAS) hits. (A) Chromosome 12q13.2 regional association plot. (B) Chromosome 17q21.1 regional association plot. The top panel highlights all genes in the region. The marginally associated TWAS genes are shown in blue, and the jointly significant genes are shown in green. The bottom panel show a regional Manhattan plot of the genome-wide association study (GWAS) data before (grey) and after (blue) conditioning on the predicted expression of the green genes.
Integrating peripheral blood eQTL for external validation
Compared with the whole blood eQTL dataset from GTEx, peripheral blood eQTL yielded fewer testable genes (n=2454) (online supplemental table S4). Across all common testable genes (n=1686), the effect sizes of those genes between GTEx and NTR were highly correlated (R=0.60, p<0.001) (online supplemental figure S2). Among the 47 TWAS significant genes, there were 15 testable genes in the replication stage. Finally, TWAS analysis using the independent peripheral blood eQTL data successfully validated 11 genes as RA risk genes (online supplemental table S5).
Comparison of TWAS with different gene-based methodologies
Venn plots illustrated the number of nominal significant genes achieved by the four methodologies. The cross-tissues test achieved more significantly associated genes, showing 67.6%, 148.7% improvement compared with FUSION and SMR, respectively (online supplemental figure S3A). Consistently, such improvement was observed within FDR-corrected thresholds. Among the 47 identified TWAS-significant associations, 85.1% (40/47) were recovered by SMR or MAGMA after FDR correction (online supplemental figure S3B).
The top 20 most significant gene set enrichment results based on (A) TWAS-based gene set enrichment analysis (TWAS-GSEA) (B) MAGMA-based gene set enrichment analysis. nine gene-sets shared between TWAS-GSEA and MAGMA-based gene set enrichment analysis with FDR <0.05 were marked in red. FDR, false discovery rate; MAGMA, multimarker analysis of genomic annotation; TWAS, transcriptome-wide association study.
Differentiation analysis in RA and controls
To explore whether these 47 risk genes identified in TWAS analysis are differentially expressed between RA patients and health controls, we evaluated their mRNA expression levels in whole blood from RA patients and healthy controls. In lines with TWAS results, eight genes including PTPN22, PUS10, PPIP5K2, ZNF322, IRF5, PSMD5, YAF2 and TMEM50B were remarkably overexpressed, while seven genes including PGBD1, PHF19, PGAP3, GSDMB, ORMDL3, TYK2 and UBASH3A were decreased in RA patients compared with healthy controls (online supplemental table S6).
Pathway enrichment analysis
Across all the 828 gene-sets with at least 10 cis-heritable genes, 10 gene-sets were significantly enriched after multiple testing correction (online supplemental table S7). MAGMA-based gene-set analysis suggested that 45 gene-sets were significantly enriched (PFDR<0.05) (online supplemental table S8). Among the top 20 enriched pathways, nine significant gene-sets were consistent in both TWAS-GSEA and MAGMA (figure 3). The results provided strong evidence of shared aetiology among autoimmune thyroid disease, asthma, type I diabetes and RA and pointed out several crucial immune-related signalling such as the Th1/Th2, intestinal immune network for IGA production, antigen processing and presentation in the pathogenesis of RA.
Discussion
With the currently released largest RA GWAS dataset, we systematically estimated the associations between genetically predicted gene expressions and RA risk. A total of 47 susceptibility genes for RA were prioritised through a two-stage TWAS design, and 13 of which resided in the novel loci. Together with in silico analysis deepened our understanding of genetic architecture and comorbidity aetiology of RA.
According to our results, the cross-tissue TWAS approach can effectively achieve more significantly associated genes compared with other two single-tissue gene-based methodologies. A recent cross-tissue TWAS has discovered two novel carcinogenic susceptibility genes for lung cancer, of which the design is similar to that of our study.26 Liu et al27 conducted a cross-tissue TWAS for pancreatic cancer, in which they revealed 13 significant gene-level associations at an FDR below 0.05, including six new susceptibility genes. Focused on RA, there has been a study of TWAS that generates genetic expression matrix in four tissues from GTEx dataset separately with FUSION software.28 The authors reported a total of 692 significant TWAS genes with p value less than 0.05. To circumvent the drawbacks of single-tissue analysis, we adopted cross-tissue analyses for discovery strategy in order to identify more reliable genes. Additionally, with expanded data from GTEx project and a more strict threshold of significance, our study is likely to achieve stable and accurate results.
Several TWAS-significant genes with strong evidence from previous functional studies were located at the known RA loci. For example, hyperactivity of PTPN22 in vitro might lead to produce reactive oxygen species, driving RA through abnormal inflammatory response and joint damage.29 PTPN22 has also been suggested as a preclinical molecular signature for RA.30 Existing evidence has proved that a part of novel RA genes exhibits immune-related features. For instance, ZNF322 at 6p22.2, which belongs to the zinc finger protein family, has been reported to act as a transcriptional activator in MAPK signalling pathways.31 Apart from this, ZNF322 can regulate the expression of two cell-cycle genes (P27 and CDK2), which were characterised as crucial regulators involved in synoviocytes proliferation.32–35 SPNS1, which encodes Spinster homologue 1 protein, is a transmembrane protein that can modulate the metabolism of autophagic lysosomal and critically associate with cellular ageing and survival.36 37 Notably, SPNS1 has recently been considered as a core component of T cell receptor signalling.38 A slew of studies have shown that TRIM10 and TRIM27 were involved in several inflammation processes.39–41 TRIM10 exhibited a lower expression in patients with systemic lupus erythematosus than in healthy individuals, which negatively induced the IFN/JAK/STAT signalling pathway through impacting the interaction between IFNAR1 and TYK2.39 Of note is that TRIM10 has been reported to be associated with systemic juvenile idiopathic arthritis.42 Moreover, Liu et al43 demonstrated that TRIM27 modulated the proliferation of mesangial cell in kidneys of lupus nephritis. It has been suggested that knockdown of TRIM27 inhibited the endothelial cells injuries in lupus nephritis via the FoxO1 signalling pathway.44
Numerous observational studies have observed the coexistence of RA and other immune-mediated diseases, consistent with significant genetic correlations between them.45–48 In our work, both TWAS-based and MAGMA-based enrichment analysis also illustrated the shared genetic determinants between autoimmune thyroid disease, asthma, type I diabetes mellitus and RA. Speculatively, there are common immune regulatory mechanisms among these diseases. Insights gained from these results will provide further therapeutic directions for immune-mediated diseases.
Although TWAS has advantages in improving statistical power and avoiding reverse causality, some caveats should be noted in this study. First, not all genes can be captured due to the criterion of significant cis-heritability genes in the TWAS analyses, and those SNPs influencing RA but are independent of cis expression will be ignored. Second, the genetically predicted models were measured in multiple tissues but not in the biologically relevant tissues such as synovial tissues and immune cells. As high-throughput data continue to be released for more cell types and tissues, as well as a larger RA GWAS dataset of 276 020 samples from five ancestral groups,49 cross-tissue association analysis promises to show even better effectiveness and provides greater insights for RA genetics. Finally, the underlying mechanisms of association signals have not been validated through experimental techniques.
In summary, we reveals 13 novel susceptibility genes whose genetically predicted expression is associated with risk of RA, providing new insights into the underlying genetic architecture of RA. However, further functional studies are still needed to elucidate the underlying biological activities of these significant signals.
Data availability statement
The datasets used and/or analysed during the current study are available from the corresponding authors on reasonable request.
Ethics statements
Patient consent for publication
Ethics approval
Not applicable.
Acknowledgments
We would like to thank the computing platform provided by Medical Data Processing Center of the School of Public Health in Anhui Medical University and all the study participants and research staff for their contributions and commitment to the present study.
References
Supplementary materials
Supplementary Data
This web only file has been produced by the BMJ Publishing Group from an electronic file supplied by the author(s) and has not been edited for content.
Supplementary Data
This web only file has been produced by the BMJ Publishing Group from an electronic file supplied by the author(s) and has not been edited for content.
Footnotes
JN, PW and K-JY are joint first authors.
Contributors H-FP coordinated and supervised the overall study. JN and PW conceived, designed the study and drafted the manuscript. K-JY and JN performed data imputation, the association analysis and bioinformatic analysis. X-KY, HC and CS commented and G-CW helped in revising the manuscript. All authors reviewed and approved the manuscript. H-FP and JN is responsible for the overall content as the guarantor.
Funding This study was funded by grants from the Natural Science Foundation of Anhui Province (2108085QH361, 2108085Y26), National Natural Science Foundation of China (81872687) and Natural Science Foundation of Anhui Medical University (2020xkj011) and Research Fund of Anhui Institute of Translational Medicine (2021zhyx-B04, 2021zhyx-C28).
Competing interests None declared.
Provenance and peer review Not commissioned; externally peer reviewed.
Supplemental material This content has been supplied by the author(s). It has not been vetted by BMJ Publishing Group Limited (BMJ) and may not have been peer-reviewed. Any opinions or recommendations discussed are solely those of the author(s) and are not endorsed by BMJ. BMJ disclaims all liability and responsibility arising from any reliance placed on the content. Where the content includes any translated material, BMJ does not warrant the accuracy and reliability of the translations (including but not limited to local regulations, clinical guidelines, terminology, drug names and drug dosages), and is not responsible for any error and/or omissions arising from translation and adaptation or otherwise.