# Investigating the GWAS-Implicated Loci for Rheumatoid Arthritis in the Pakistani Population

Authors: Muhammad Muaaz Aslam; Peter John; Kang-Hsien Fan; Attya Bhatti; Wajahat Aziz; Bashir Ahmed; Eleanor Feingold; F. Yesim Demirci; M. Ilyas Kamboh
Journal: Disease Markers
 (2020)

---

## Abstract

Rheumatoid arthritis (RA) is a complex and multifactorial autoimmune disorder with the involvement of multiple genetic and environmental factors. Genome-wide association studies (GWAS) have identified more than 50 RA genetic loci in European populations. Given the anticipated overlap of RA-relevant genes and pathways across different ethnic groups, we sought to replicate 58 GWAS-implicated SNPs reported in Europeans in Pakistani subjects. 1,959 unrelated subjects comprising 1,222 RA cases and 737 controls were collected from three rheumatology facilities in Pakistan. Genotyping was performed using iPLEX or TaqMan® methods. A total of 50 SNPs were included in the final association analysis after excluding those that failed assay design/run or postrun QC analysis. Fourteen SNPs (LINC00824/rs1516971, PADI4/rs2240336, CEP57/rs4409785, CTLA4/rs3087243, STAT4/rs13426947, HLA-B/MICA/rs2596565, C5orf30/rs26232, CCL21/rs951005, GATA3/rs2275806, VPS37C/rs595158, HLA-DRB1/rs660895, EOMES/rs3806624, SPRED2/rs934734, and RUNX1/rs9979383) were replicated in our Pakistani sample at false discovery rate (FDR) of <0.20 with nominal p values ranging from 4.73E-06 to 3.48E-02. Our results indicate that several RA susceptibility loci are shared between Pakistani and European populations, supporting the role of common genes/pathways.


---

## Body

## 1. Introduction

Rheumatoid arthritis (RA) is a complex autoimmune disorder due to the involvement of many genetic, environmental, and stochastic elements in its etiology [1, 2]. It is characterized by continuous inflammation associated with altered expression of different proinflammatory (TNFα, IL-1, IL-17) and anti-inflammatory (IL-4, IL-10, IL-13, IL-35) cytokines and activation of B and T cells, which leads to the destruction of synovial joints and ultimately physical disability [3, 4]. Worldwide, RA affects 0.5 to 1% of the population [5]. RA prevalence does not differ significantly between rural and urban areas [6]. It can affect both sexes at any age, but data suggest that women are three times more prone than men [7]. There is a significant genetic susceptibility to RA with 50-60% heritability estimates [8]. Twin studies also revealed high concordance rates in monozygotic twins (12.3-15.4%) relative to dizygotic twins (3.5%), strongly indicating the presence of a genetic component in RA [9].During the last two decades, many genetic studies, including genome-wide association studies (GWAS), have been conducted to understand the genetic basis of RA and this has led to the identification of more than 50 risk loci [10]. Besides HLA-DRB1, many non-HLA loci have been implicated that are involved in multiple functions/pathways including immune activation (NF-κB signaling, JAK-STAT pathway, regulation, and activation of CD4+T-cells, complement activation), fibroblast differentiation and dedifferentiation, and bone modeling and repair [11, 12].Thus far, the major focus of genetic studies on RA has largely been on Europeans or European-derived populations, with only limited data available on other populations [13], especially in South Asians [14]. In an effort to expand genetic studies on RA in South Asians and to further delineate its genetic basis, we sought to replicate 58 genome-wide significant single nucleotide polymorphisms (SNPs) from 58 loci previously implicated in Europeans or European-derived populations [15–18] in Pakistani population where there is the paucity of such data.


## 2. Material and Methods

### 2.1. Study Subjects

A total of 1,959 subjects (1,222 cases, 737 controls) were included in our study. Blood samples and relevant clinical information were collected from three major public or private rheumatology clinics in Pakistan: Pakistan Institute of medical sciences (PIMS), Military hospital, and Rehmat Noor Clinic. All cases included in this study (meanage±SD=43.1±12.33, 78.6% women) were clinically diagnosed by rheumatologists and met the American College of Rheumatology (ACR) 1987 classification criteria for RA [19]. All controls included in this study (meanage±SD=40.7±12.49, 39.5% women) were recruited from the general population and were free of any autoimmune disease at the time of recruitment. A screening questionnaire was filled out and a written informed consent was obtained from each subject at the time of the recruitment. All blood samples were collected in EDTA tubes to avoid coagulation and processed shortly after the collection. The study was approved by the Institutional review board (IRB) of the University of Pittsburgh, USA (IRB no. PRO12110472).


### 2.2. Genomic DNA Extraction

Genomic DNA was extracted from whole blood using either the standard phenol-chloroform extraction method or GeneJET Whole Blood Genomic DNA Purification (Thermo Scientific USA) and quantified using the NanoDrop™ 2000 spectrophotometer (Thermo Scientific USA).


### 2.3. Genotyping

A total of 58 RA-associated genome-wide significant SNPs (p<5E−08) was selected from the previously published GWAS in subjects of European ancestry (Table 1). SNP genotyping was performed using either the TaqMan® (Applied Biosystems, Thermo Fisher Scientific) or iPLEX® Gold (Agena Bioscience) methods and following the manufacturer’s design/order instructions and protocols. After the thermal cycling of the TaqMan® assays and DNAs on 384-well plates, the endpoint fluorescence reading was performed on a QuantStudio™ 12K Flex system (Applied Biosystems, Thermo Fisher Scientific). The iPLEX® Gold genotyping was performed in the Core laboratories of the University of Pittsburgh. 18% replicates were used to test genotyping consistency.Table 1
List of selected GWAS-implicated RA SNPs examined in this study.
SNPLocusChromosomal location (GRCh38)Variant typeGeneReferencers22281451q21154454494Missense variantIL6R[15]rs21053251q25.1173380586Intron variantLOC100506023[16]rs28434011p362596694Intron variantMMEL1[15]rs284113521p34.3378129073 prime UTR variantMTF1[16]rs8832201p3438151199Intron variantLOC105378654[15]rs24766011p13113834946Missense variantPTPN22[15]rs67325652q13110850255Intron variantACOXL[16]rs116769222q11100190478IntergenicAFF3[17]rs67152842q33.1201289674Intron variantALS2CR12[16]rs30872432q33203874196Downstream variantCTLA4[17]rs346959442p1660897715Intron variantREL[15]rs9347342p1465368452Intron variantSPRED2[17]rs134269472q32191068528Intron variantSTAT4[15]rs38066243p24.127723132Upstream variantEOMES[16]rs98268283q22.3136683218Intron variantSTAG1[16]rs44523133p24.317005540Intron variantPLCL2[16]rs133155913p1458571114Intron variantFAM107A[17]rs8740404p1526106575Upstream variantRBPJ[17]rs26640354p1148218822Intron variantTEC[16]rs716241195q1156144903Intron variantANKRD55[15]rs262325q21103261019Intron variantC5orf30[17]rs30930236q27167120802Intron variantCCR6[16]rs22340676p21.3136387877Upstream variantETV7[16]rs69202206q23137685367Upstream variantTNFAIP3[15]rs25965656p21.3331385552Upstream variantHLA-B/MICA[18]rs69100716p2132315077Intron variantTSBP1[17]rs6608956p2132609603Intron variantHLA-DRB1[14]rs42727q21.2926075153 prime UTR variantCDK6[16]rs104886317q32128954129Downstream variantIRF5[17]rs672504507p15.128135367Intron variantJAZF1[16]rs6783478q22.3101451374Upstream variantGRHL2[16]rs15169718q24.21128529854Intron variantLINC00824[16]rs9987318q21.1380183160Intron variantTPD52[16]rs9510059p1334743684IntergenicCCL21[17]rs1276437810q2162040245Intron variantARID5B[15]rs227580610p148053377Intron variantGATA3[15]rs70677810p156056986Intron variantIL2RA[17]rs59515811q1261142109Intron variantVPS37C[15]rs440978511q2195578258Intergenic/unknownCEP57[16]rs77312512q13.256001170Intron variantSUOX[16]rs195089714q24.168293424Intron variantRAD51B[16]rs804308515q1438535939Intron variantRASGRP1[15]rs802689815q2369699078Intergenic/unknownTLE3[15]rs1333017616q2485985481Intergenic/unknownIRF8[15]rs478040116p13.1311745470Upstream variantTXNDC11[16]rs1293640917q1239887396Intergenic/unknownIKZF3[15]rs3453644319p1310352442Missense variantTYK2[15]rs481048520q1346119308Intron variantCD40[17]rs997938321q2235343463Intron variantRUNX1[15]rs189359221q22.342434957Intron variantUBASH3A[16]rs90968522q13.139351666Intron variantSYNGR1[16]rs22403361p3617347907Intron variantPADI4[16]rs101757982p23.130226728Upstream variantLBH[16]rs96856711q12.261828092Intron variantFADS2[16]rs1077462412q24.12111395984Upstream variantLINC02356[16]rs960361613q14.1139793932Downstream variantCOG6[16]rs7319405821q22.1133391982Intergenic/unknownIFNGR2[16]rs283451221q2234539301Intron variantRCAN1[15]


### 2.4. Statistical Analysis

Concordance to Hardy-Weinberg Equilibrium (HWE) was tested using the Chi-square goodness of fit test. Departure from HWE was considered at arbitrarilyp<1E−05. Logistic regression using an additive model and minor allele as the effect allele was employed for case-control association analysis using sex and age as covariates. The Benjamin Hochberg false discovery rate (FDR) was applied to correct for multiple testing [20]. p<0.05 was considered as suggestive evidence of association and FDR (q value) of <0.20 as statistically significant as used in previous reports [21, 22]. All analyses were implemented in R, version 3.4.4.


### 2.5. Functional Annotations

To evaluate the potential biological significance of reported genome-wide significant SNPs, we used the Genotype-Tissue Expression (GTEx) database (https://gtexportal.org/home/) to search for expression quantitative trait loci (eQTL) in RA-relevant tissues and whole blood. We also used the RegulomeDB online database (http://regulome.stanford.edu/) to determine possible regulatory functions of the SNPs located in noncoding regions.


## 2.1. Study Subjects

A total of 1,959 subjects (1,222 cases, 737 controls) were included in our study. Blood samples and relevant clinical information were collected from three major public or private rheumatology clinics in Pakistan: Pakistan Institute of medical sciences (PIMS), Military hospital, and Rehmat Noor Clinic. All cases included in this study (meanage±SD=43.1±12.33, 78.6% women) were clinically diagnosed by rheumatologists and met the American College of Rheumatology (ACR) 1987 classification criteria for RA [19]. All controls included in this study (meanage±SD=40.7±12.49, 39.5% women) were recruited from the general population and were free of any autoimmune disease at the time of recruitment. A screening questionnaire was filled out and a written informed consent was obtained from each subject at the time of the recruitment. All blood samples were collected in EDTA tubes to avoid coagulation and processed shortly after the collection. The study was approved by the Institutional review board (IRB) of the University of Pittsburgh, USA (IRB no. PRO12110472).


## 2.2. Genomic DNA Extraction

Genomic DNA was extracted from whole blood using either the standard phenol-chloroform extraction method or GeneJET Whole Blood Genomic DNA Purification (Thermo Scientific USA) and quantified using the NanoDrop™ 2000 spectrophotometer (Thermo Scientific USA).


## 2.3. Genotyping

A total of 58 RA-associated genome-wide significant SNPs (p<5E−08) was selected from the previously published GWAS in subjects of European ancestry (Table 1). SNP genotyping was performed using either the TaqMan® (Applied Biosystems, Thermo Fisher Scientific) or iPLEX® Gold (Agena Bioscience) methods and following the manufacturer’s design/order instructions and protocols. After the thermal cycling of the TaqMan® assays and DNAs on 384-well plates, the endpoint fluorescence reading was performed on a QuantStudio™ 12K Flex system (Applied Biosystems, Thermo Fisher Scientific). The iPLEX® Gold genotyping was performed in the Core laboratories of the University of Pittsburgh. 18% replicates were used to test genotyping consistency.Table 1
List of selected GWAS-implicated RA SNPs examined in this study.
SNPLocusChromosomal location (GRCh38)Variant typeGeneReferencers22281451q21154454494Missense variantIL6R[15]rs21053251q25.1173380586Intron variantLOC100506023[16]rs28434011p362596694Intron variantMMEL1[15]rs284113521p34.3378129073 prime UTR variantMTF1[16]rs8832201p3438151199Intron variantLOC105378654[15]rs24766011p13113834946Missense variantPTPN22[15]rs67325652q13110850255Intron variantACOXL[16]rs116769222q11100190478IntergenicAFF3[17]rs67152842q33.1201289674Intron variantALS2CR12[16]rs30872432q33203874196Downstream variantCTLA4[17]rs346959442p1660897715Intron variantREL[15]rs9347342p1465368452Intron variantSPRED2[17]rs134269472q32191068528Intron variantSTAT4[15]rs38066243p24.127723132Upstream variantEOMES[16]rs98268283q22.3136683218Intron variantSTAG1[16]rs44523133p24.317005540Intron variantPLCL2[16]rs133155913p1458571114Intron variantFAM107A[17]rs8740404p1526106575Upstream variantRBPJ[17]rs26640354p1148218822Intron variantTEC[16]rs716241195q1156144903Intron variantANKRD55[15]rs262325q21103261019Intron variantC5orf30[17]rs30930236q27167120802Intron variantCCR6[16]rs22340676p21.3136387877Upstream variantETV7[16]rs69202206q23137685367Upstream variantTNFAIP3[15]rs25965656p21.3331385552Upstream variantHLA-B/MICA[18]rs69100716p2132315077Intron variantTSBP1[17]rs6608956p2132609603Intron variantHLA-DRB1[14]rs42727q21.2926075153 prime UTR variantCDK6[16]rs104886317q32128954129Downstream variantIRF5[17]rs672504507p15.128135367Intron variantJAZF1[16]rs6783478q22.3101451374Upstream variantGRHL2[16]rs15169718q24.21128529854Intron variantLINC00824[16]rs9987318q21.1380183160Intron variantTPD52[16]rs9510059p1334743684IntergenicCCL21[17]rs1276437810q2162040245Intron variantARID5B[15]rs227580610p148053377Intron variantGATA3[15]rs70677810p156056986Intron variantIL2RA[17]rs59515811q1261142109Intron variantVPS37C[15]rs440978511q2195578258Intergenic/unknownCEP57[16]rs77312512q13.256001170Intron variantSUOX[16]rs195089714q24.168293424Intron variantRAD51B[16]rs804308515q1438535939Intron variantRASGRP1[15]rs802689815q2369699078Intergenic/unknownTLE3[15]rs1333017616q2485985481Intergenic/unknownIRF8[15]rs478040116p13.1311745470Upstream variantTXNDC11[16]rs1293640917q1239887396Intergenic/unknownIKZF3[15]rs3453644319p1310352442Missense variantTYK2[15]rs481048520q1346119308Intron variantCD40[17]rs997938321q2235343463Intron variantRUNX1[15]rs189359221q22.342434957Intron variantUBASH3A[16]rs90968522q13.139351666Intron variantSYNGR1[16]rs22403361p3617347907Intron variantPADI4[16]rs101757982p23.130226728Upstream variantLBH[16]rs96856711q12.261828092Intron variantFADS2[16]rs1077462412q24.12111395984Upstream variantLINC02356[16]rs960361613q14.1139793932Downstream variantCOG6[16]rs7319405821q22.1133391982Intergenic/unknownIFNGR2[16]rs283451221q2234539301Intron variantRCAN1[15]


## 2.4. Statistical Analysis

Concordance to Hardy-Weinberg Equilibrium (HWE) was tested using the Chi-square goodness of fit test. Departure from HWE was considered at arbitrarilyp<1E−05. Logistic regression using an additive model and minor allele as the effect allele was employed for case-control association analysis using sex and age as covariates. The Benjamin Hochberg false discovery rate (FDR) was applied to correct for multiple testing [20]. p<0.05 was considered as suggestive evidence of association and FDR (q value) of <0.20 as statistically significant as used in previous reports [21, 22]. All analyses were implemented in R, version 3.4.4.


## 2.5. Functional Annotations

To evaluate the potential biological significance of reported genome-wide significant SNPs, we used the Genotype-Tissue Expression (GTEx) database (https://gtexportal.org/home/) to search for expression quantitative trait loci (eQTL) in RA-relevant tissues and whole blood. We also used the RegulomeDB online database (http://regulome.stanford.edu/) to determine possible regulatory functions of the SNPs located in noncoding regions.


## 3. Results

A total of 1,222 unrelated RA cases and 737 controls were recruited for this research study. The prevalence of RA was higher in females (78%) than males (22%), supporting the earlier data that females are more prone to RA [14].Eight of the 58 genotyped SNPs failed the QC (quality control) during assay design (either iPLEX® Gold/TaqMan® or both). The genotype distribution of all QC-passed 50 SNPs adhered to the HWE. The association analyses results in our Pakistani sample are presented in Table2. Fourteen SNPs showed nominal significance at p<0.05 and FDR of ≤0.20.Table 2
Association analysis results for GWAS-implicated RA SNPs in the Pakistani population.
Gene/SNPMajor alleleMinor alleleMAFReported GWASp valueNominalp value in PakistanisOR (95% CI)FDR (q)LINC00824/rs1516971TC0.0993.20E-114.73E-060.57 (0.45, 0.73)2.37E-04PADI4/rs2240336CT0.4525.9E−95.00E-050.74 (0.64, 0.86)1.25E-03CEP57/rs4409785TC0.2623.60E-091.03E-031.33 (1.12, 1.58)1.54E-02CTLA4/rs3087243AG0.4221.2E−81.23E-031.28 (1.1, 1.49)1.54E-02STAT4/rs13426947GA0.2387.2E−102.59E-031.31 (1.1, 1.56)2.59E-02HLA-B/MICA/rs2596565GA0.1199.26E-094.53E-031.4 (1.11, 1.77)3.77E-02C5orf30/rs26232CT0.1814.10E-083.73E-020.83 (0.69, 0.99)1.52E-01CCL21/rs951005AG0.2483.90E-102.75E-020.82 (0.69, 0.98)1.52E-01GATA3/rs2275806AG0.3974.6E−84.02E-021.17 (1.01, 1.36)1.52E-01VPS37C/rs595158AC0.3203.4E−83.10E-021.19 (1.02, 1.38)1.52E-01HLA-DRB1/rs660895AG0.111<1E-3004.02E-021.27 (1.01, 1.59)1.52E-01EOMES/rs3806624GA0.2622.80E-083.54E-020.84 (0.71, 0.99)1.52E-01SPRED2/rs934734AG0.4065.30E-104.24E-021.17 (1.01, 1.35)1.52E-01RUNX1/rs9979383TC0.3005.0E−103.48E-020.84 (0.72, 0.99)1.52E-01IL2RA/rs706778TC0.4581.40E-116.21E-020.87 (0.75, 1.01)2.07E-01LOC105378654/rs883220CA0.2122.1E−81.08E-010.87 (0.73, 1.03)3.17E-01IRF5/rs10488631TC0.1964.20E-111.07E-011.16 (0.97, 1.39)3.17E-01RCAN1/rs2834512GA0.1602.1E−81.39E-010.86 (0.71, 1.05)3.86E-01ARID5B/rs12764378GA0.2414.5E−101.64E-011.13 (0.95, 1.34)4.33E-01TYK2/rs34536443GC0.0102.3E−141.73E-010.62 (0.31, 1.23)4.33E-01IL6R/rs2228145AC0.3231.3E−81.89E-010.9 (0.77, 1.05)4.50E-01PTPN22/rs2476601GA0.0167.5E−772.21E-011.47 (0.79, 2.73)4.93E-01SYNGR1/rs909685TA0.4806.40E-122.27E-011.09 (0.95, 1.26)4.93E-01MTF1/rs28411352CT0.2225.90E-092.42E-011.11 (0.93, 1.32)5.05E-01CCR6/rs3093023GA0.4251.50E-112.70E-011.09 (0.94, 1.26)5.19E-01ALS2CR12/rs6715284CG0.1852.50E-092.61E-011.11 (0.92, 1.35)5.19E-01UBASH3A/rs1893592AC0.2419.80E-093.11E-010.92 (0.77, 1.08)5.75E-01IRF8/rs13330176TA0.2434.0E−84.04E-010.93 (0.78, 1.11)6.13E-01PLCL2/rs4452313AT0.4925.2E-113.76E-011.07 (0.92, 1.24)6.13E-01ACOXL/rs6732565GA0.4449.40E-094.15E-010.94 (0.81, 1.09)6.13E-01GRHL2/rs678347AG0.4257.30E-104.17E-011.06 (0.92, 1.23)6.13E-01TLE3/rs8026898GA0.3039.16.E-143.50E-010.93 (0.79, 1.09)6.13E-01RASGRP1/rs8043085GT0.2801.4E−103.64E-011.08 (0.92, 1.27)6.13E-01REL/rs34695944TC0.1121.4E−104.05E-011.1 (0.88, 1.39)6.13E-01COG6/rs9603616CT0.2602.80E-114.47E-010.94 (0.8, 1.11)6.38E-01TPD52/rs998731TC0.4346.60E-095.04E-010.95 (0.82, 1.1)7.00E-01TNFAIP3/rs6920220GA0.1322.3E−135.33E-011.07 (0.86, 1.34)7.20E-01CD40/rs4810485GT0.2502.8E−96.03E-010.96 (0.81, 1.13)7.94E-01RAD51B/rs1950897TC0.2085.00E-087.02E-011.04 (0.86, 1.24)8.98E-01FAM107A/rs13315591TC0.0434.60E-087.25E-011.06 (0.75, 1.51)8.98E-01MMEL1/rs2843401CT0.4616.6E−97.80E-011.02 (0.88, 1.18)8.98E-01RBPJ/rs874040GC0.1591.00E-168.20E-010.98 (0.81, 1.18)8.98E-01IKZF3/rs12936409CT0.3922.8E−98.19E-010.98 (0.85, 1.14)8.98E-01LBH/rs10175798AG0.4744.20E-088.46E-011.01 (0.88, 1.17)8.98E-01LINC02356/rs10774624AG0.1386.90E-098.38E-010.98 (0.79, 1.2)8.98E-01FADS2/rs968567GA0.0841.80E-088.62E-010.98 (0.75, 1.27)8.98E-01AFF3/rs11676922TA0.4851.00E-147.91E-010.98 (0.85, 1.14)8.98E-01CDK6/rs4272AG0.0831.2E-87.64E-010.96 (0.74, 1.24)8.98E-01TEC/rs2664035GA0.2563.30E-089.54E-011 (0.84, 1.17)9.70E-01TXNDC11/rs4780401TG0.4308.70E-099.70E-011 (0.86, 1.16)9.70E-01We found the same direction of association where the minor allele showed the same risk or protective effect as compared to the previous findings on the tested allele (major or minor). In our data,LINC00824/rs1516971 showed the most significant association (p=4.73E−06). The second most significant SNP was PADI4/rs2240336 (p=5.00E−05) followed by CEP57/rs4409785 (p=1.03E−03), CTLA4/rs3087243 (p=1.23E−03), STAT4/rs13426947 (p=2.59E−03), and HLA-B-MICA/rs2596565 (p=4.53E−03). Eight other SNPs showed marginal significance: C5orf30/rs26232 (p=3.73E−02), CCL21/rs951005 (p=2.75E−02), GATA3/rs2275806 (p=4.02E−02), VPS37C/rs595158 (p=3.10E−02), HLA-DRB1/rs660895 (p=4.02E−02), EOMES/rs3806624 (p=3.54E−02), SPRED2/rs934734 (p=4.24E−02), and RUNX1/rs9979383 (p=3.48E−02). Figure 1 shows the distribution of tested SNPs across the genome where the SNPs with p value < 0.05 are labeled.Figure 1
Annotated 50 tested SNPs. SNPs withp value < 0.05 are shown above the dotted line.Next, we examined the functional significance of all 50 SNPs using the GTEx and RegulomeDB databases. Table3 shows the category summaries of RegulomeDB scores, and Table 4 shows 14 SNPs that had a RegulomeDB score of ≤3, indicating strong evidence of potential regulatory role. SNPs falling in this category have more likelihood to affect the binding of transcriptional factors. Out of these fourteen SNPs, only five (HLA-B/MICA/rs2596565, HLA-DRB1/rs660895, C5orf30/rs26232, CTLA4/rs3087243, and RUNX1/rs9979383) had p<0.05 in our association results. SYNGR1/rs909685 had the top RegulomeDB score of 1b followed by FADS2/rs968567, CD40/rs4810485, CCR6/rs3093023, HLA-B/MICA/rs2596565, and HLA-DRB1/rs660895 with a score of 1f. The latter two SNPs (HLA-B/MICA/rs2596565 and HLA-DRB1/rs660895) were significant in our sample (p=0.0045 and p=0.0402) and based on RegulomeDB, both SNPs are eQTL for HLA-DQA1. While HLA-DRB1/rs660895 affects the binding of two proteins (BCLAF1 and POLR2A), no such evidence was found for HLA-B/MICA/rs2596565. SYNGR1/rs909685 is an intronic variant, which falls in the p53decamer binding motif and affects the binding of five different proteins (MAX, MYC, PAX5, TRIM28, and EBF1). FADS2/rs968567 is also an intronic variant and it affects the binding of twenty-eight proteins, including MAX, POLR2A, NFKB1, and NFIC. The binding of these proteins is also affected by the CD40/rs4810485 variant. An intronic CCR6/rs3093023 variant falls in the Oct-1 binding motif and is eQTL for CCR6.Table 3
RegulomeDB score description.
ScoreDescriptionLikely to affect binding and linked to the expression of a gene target1aeQTL + TF binding + matched TF motif + matched DNase footprint + DNase peak1beQTL + TF binding + any motif + DNase footprint + DNase peak1ceQTL + TF binding + matched TF motif + DNase peak1deQTL + TF binding + any motif + DNase peak1eeQTL + TF binding + matched TF motif1feQTL + TF binding/DNase peakLikely to affect binding2aTF binding + matched TF motif + matched DNase footprint + DNase peak2bTF binding + any motif + DNase footprint + DNase peak2cTF binding + matched TF motif + DNase peakLess likely to affect binding3aTF binding + any motif + DNase peak3bTF binding + matched TF motifMinimal binding evidence4TF binding + DNase peak5TF binding or DNase peak6Motif hit7No data availableTable 4
Details of study SNPs (RegulomeDBScore≤3) with putative regulatory functions.
SNPChr.ScoreeQTLBound proteinMotifsSYNGR1/rs909685chr221bSYNGR1MAX, MYC, PAX5, TRIM28, EBF1p53decamerFADS2/rs968567chr111fNXF1CTCF, E2F4, EP300, ETS1, FOSL2, FOXP2, GABPA, GABPB1, GATA1, HNF4A, HNF4G, MAX, MYBL2, NFIC, NFKB1, NFYB, PML, POLR2A, RAD21, REST, SIN3A, SP1, SREBF1, SREBF2, TAF1, TCF12, YY1, ZBTB7ACD40/rs4810485chr201fCD40POLR2A, NFKB1, NFIC, MEF2C, MEF2A, IRF1, IKZF1, FOXM1, BCL3, BATFCCR6/rs3093023chr61fCCR6Oct_1HLAB/MICA/rs2596565chr61fBTN3A2, HLA-A, HLA-C, HLA-DPB1, HLA-DQA1, HLA-DQB1, HLA-DRB1, HLA-G, HLA-H, VARSLHLA-DRB1/rs660895chr61fHLA-DQA2, HLA-DQA1BCLAF1, POLR2AGRHL2/rs678347chr82aYY1, CEBPB, EP300, FOXA2, JUND, STAT3, TCF12, USF2, SIN3AHNF3, Foxa2, FOXF2, Freac-2, Freac-4, FOXA1, HNF3, FOXC1PTPN22/rs2476601chr12bFOS, STAT3Gabpa, Etv1, Elk1, Erg, Ehf, PU.1LBH/rs10175798chr22bATF2, FOXM1, RUNX3, NFIC, SPI1, MEF2A, BATFSTAT1:STAT1RBPJ/rs874040chr42bGATA2, PML, TAL1MEF-2C5orf30/rs26232chr52bBHLHE40, CHD1, EP300, FOXM1, MAZ, NFATC1, NFIC, RUNX3, TBL1XR1, TBP, USF2, ATF2, MEF2A, MEF2C, NFKB1ICSBPRUNX1/rs9979383chr213aMYC, NFKB1E47, FIGLA, ID4, MESP1, SNAI2CTLA4/rs3087243chr23aMAX, FOS, STAT3BbxIRF5/rs10488631chr73aGATA2NanogGTEx data eQTL showed the lowest p-eQTL (2.59E-34) for FADS2/rs968567 that affects FADS2 gene expression in transformed fibroblast cells; the same variant also affects TMEM258 (p−eQTL=3.43E−08), ZBTB3 (p−eQTL=0.0168), MYRF (p−eQTL=0.0207), DAGLA (p−eQTL=0.0246), and RAB3IL1 (p−eQTL=0.0247). In transformed fibroblast cells, LINC00824/rs1516971 was a weak eQTL (p=0.0157) for RP11-89M16.1. No data was available for CTLA4/rs3087243 in GTEx.STAT4/rs13426947 was eQTL for STAT4, AC005540.3, and RP11-647K16.1 in EBV-transformed lymphocyte cells. HLA-B/MICA/rs2596565 was eQTL for multiple immune-related genes including C4A, C4B, HLA-S, HLA-B, HLA-C, and MICA in all three tested tissues. C5orf30/rs26232 showed the strongest eQTL (p−eQTL=1.72E−19) with PPIP5K2 only in whole blood. Additional details of GTEx data are given in supplementary Table S1.


## 4. Discussion

There have been a number of genome-wide association studies on RA over the past decade that have resulted in the identification of many RA susceptibility loci, thus improving our understanding of the complex genetic underpinning of RA [15–17]. Most of the recently published GWAS have been conducted on Europeans and North Americans, which have identified many new RA risk loci and replicated and confirmed the previously reported putative risk loci. Many of these reported risk loci are shared among different ethnic groups, while some are specific to certain populations [23]. Replication studies across different ethnic groups are necessary to better understand the globe-wide population specificity of these established RA risk loci and guide future directions. Since South Asians, including the Pakistani population, have significant European ancestry [24, 25], we chose European originated GWAS-implicated SNPs for replication in Pakistanis. For this purpose, we enrolled 1,959 unrelated RA cases and controls from two public and private Rheumatology facilities, where people from different demographic regions of Pakistan visit for treatment, and we successfully genotyped 50 GWAS—implicated SNPs in those subjects.In our data set, rs1516971 was the most significant SNP (p=4.73E−06), which is an intronic variant in the LINC00824 gene at chromosome 8q24.21. Previously, this variant has also been reported in a trans-ethnic association study of RA (p=1.3E−10) [16]. Another SNP (rs6651252) from the same region, which is in complete LD (r2=1.00) with rs1516971, has also been reported to be associated with RA and Crohn’s disease [15, 26, 27]. GTEx expression data suggests that rs1516971 affects the expression of RP11-89M16.1 at chromosome 8 in transformed fibroblast cells (p−eQTL=0.0157).The second most significant SNP in our data was rs2240336 (p=5.00E−05), which is present in the 9th intron of PADI4 at chromosome 1p36. PADI4 is one of the most important RA susceptibility loci in multiple ethnic groups, including Europeans, Asians, and Latin Americans [28, 29]. Another SNP in this region, rs2301888 that is in strong LD (r2=0.8) with rs2240336, has also been reported to be associated with RA in Europeans and Koreans [30]. However, a study on Iranian population found no significant association of two other SNPs (rs11203367 and rs874881) in the PADI4 region with RA, which are not in strong LD (r2<0.8) with rs2240336 [31]. GTEx expression data showed that rs2240336 is eQTL for RP4-798A10.2 (p−eQTL=0.0191) and PADI3 (p−eQTL=0.0203) in EBV transformed lymphocytes. The rs2240336 SNP also controls the expression of FBXO42 (p−eQTL=0.00526) in transformed fibroblasts. The third most important SNP in our data is CEP57/rs4409785 (p=1.03E−03), which is an intergenic variant on chromosome 11q21. This SNP has also been reported to be associated with other immune-related diseases such as multiple sclerosis, vitiligo, and Graves’ disease [32–34]. In GTEx expression data, rs4409785 was eQTL for JRKL (p−eQTL=0.037) in EBV-transformed lymphocytes, and for AP001877.1 (p−eQTL=0.0188), and MAML2 (p−eQTL=0.034) in transformed fibroblast cells.Among the reported genome-wide significant SNPs that did not show significant association in our sample, the top functional ones based on the GTEx and RegulomeDB database includedSYNGR1/rs909685, FADS2/rs968567, CD40/rs4810485, CCR6/rs3093023, and IKZF3/rs12936409. SYNGR1/rs909685 has a RegulomeDB score of 1b and is a strong eQTL for SYNGR1 in EBV-transformed lymphocytes (p−eQTL=3.74E−14) and whole blood (p−eQTL=5.68E−21). FADS2/rs968567 has a RegulomeDB score of 1f and is eQTL for FADS2 in EBV-transformed lymphocytes (p−eQTL=2.05E−05), whole blood and (p−eQTL=1.86E−71) transformed fibroblasts (p−eQTL=2.59E−34). CD40/rs4810485 also has a high RegulomeDB score of 1f and is strong eQTL for CD40 in whole blood (p−eQTL=1.33E−09) and transformed fibroblasts (p−eQTL=8.58E−11) but weak eQTL for CD40 in EBV-transformed lymphocytes (p−eQTL=0.0272). CCR6/rs3093023 with RegulomeDB score of 1f was the strongest eQTL for RNASET2 in transformed fibroblasts (p−eQTL=1.08E−16) and whole blood (p−eQTL=1.31E−15) with no data reported for EBV-transformed lymphocytes in GTEx. IKZF3/rs12936409 was a strong eQTL for ORMDL3 in all three tested tissues. In GTEx expression data, rs12936409 was found to control the expression of GSDMA, GSDMB, and ORMDL3 in all tested cells and whole blood.


## 5. Conclusions

We were able to successfully replicate 14 of 50 SNPs selected from previously published GWAS results in European populations in our unique Pakistani population. These findings suggest that there is a sharing of RA risk loci among different population groups. A weakness of our study is the availability of relatively small-sized sample, which may have prevented 36 SNPs from achieving the significance threshold, although they showed a trend for the same directional effects as the reported ones. Further studies using larger samples may help to identify and replicate more RA risk loci in the Pakistani population.



---

*Source: 1910215-2020-08-03.xml*