Missense variant in interleukin-6 signal transducer identified as susceptibility locus for rheumatoid arthritis in Chinese patients.

Objectives: This study aims to uncover variants of large effect size and allele frequency below 5% by sequencing all extant genes associated with rheumatoid arthritis (RA) in a homogeneous patient cohort. Patients and methods: This retrospective study was conducted between January 2001 and December 2017. We selected Chinese RA patients positive for anti-citrullinated peptide antibody (ACPA). All the 128 known candidate genes identified through genome-wide association studies were sequenced in 48 RA patients (15 males, 33 females; mean age 53.32±8.98 years; range, 32 to 75 years) and 45 controls (11 males, 34 females; mean age 32.18±9.54; range, 21 to 57 years). The exonic regions of these genes were sequenced. The resultant data were analyzed for association using single variant association and pathway-based association enrichment tests. The genetic burden due to low-frequency variants was assessed with the C-alpha test. The candidate variants that showed significant association were validated in a larger cohort of 500 RA cases (71 males, 429 females; mean age 48.6±12.2 years; range, 24 to 92 years) and 500 controls (66 males, 434 females; mean age 32.3±10.1 years; range, 21 to 73 years).
Results: Thirty-nine variants in 21 genes were identified using single variant association analysis and C-alpha test, with stepwise filtering. Among these, the missense variant in interleukin-6 signal transducer (IL-6ST) 5:55260065 (p.Cys47Phe) was significantly associated with RA in Chinese patients in Singapore.
Conclusion: Our results suggest that a mutation in IL-6ST (5:55260065) confers risk of RA in Chinese patients in Singapore.

Introduction

Rheumatoid arthritis (RA) is a chronic systemic inflammatory disease manifesting as joint pain, swelling, and destruction. The discovery of the association of RA and human leukocyte antigen (HLA) Dw4 suggested a genetic basis for the pathogenesis.[1] Subsequently, most of the disease-susceptibility genetic variants have been identified through genome-wide association studies (GWASs). Hitherto, more than 100 genes have been implicated in the pathogenesis of RA.[2]

The heritability of RA is estimated to be 55%,[3] but only 65% of this, or 36% of the total disease liability, can be explained by summing the effects of the single nucleotide polymorphisms (SNPs).[4] This missing heritability has been explained by rare variants of large effect size.[5] However, sequencing the exons of 25 genes in six autoimmune diseases failed to reveal any disease-associated rare variants.[6] Diogo et al.[7] sequenced the exons of 25 RA-associated genes in 500 RA cases and 650 controls of European ancestry and found two susceptibility variants, interleukin (IL)-2 receptor alpha and IL-2 receptor beta. Bang et al.[8] failed to identify any rare variants after sequencing the exons of 398 genes, including 106 known RA loci, in 1,217 Korean RA patients and 717 controls. Li et al.[9] sequenced the exomes of 58 RA patients and 66 controls and identified five susceptibility genes, NCR3LG1, RAP1GAP, CHCHD5, HIPK2, and DIAPH2.

We decided to sequence the exons of all 128 candidate genes known to be associated with RA together with their flanking non-coding regions, in a homogeneous patient cohort. We were interested in identifying novel risk variants of large effect size. We obtained these genes for targeted sequencing from a meta-analysis of 22 studies[10] and the RAvariome database.[11] Thus, in this study, we aimed to uncover variants of large effect size and allele frequency below 5% by sequencing all extant genes associated with RA in a homogeneous patient cohort.

Patients and Methods

This retrospective study was conducted at Tan Tock Seng Hospital between January 2001 and December 2017. The RA patients and healthy controls were derived from our prospective RA Registry and Healthy Control Tissue Bank, respectively. All participants were at least 21 years of age at study entry and fulfilled the 1987 American College of Rheumatology or the American College of Rheumatology/ European League Against Rheumatism 2010 classification criteria for RA.[12,13] We selected patients of Chinese descent who were positive for the anti-cyclic citrullinated peptide antibody (ACPA) to minimize heterogeneity.[14] We sequenced the 128 genes in 48 RA patients (15 males, 33 females; mean age 53.32±8.98 years; range, 32 to 75 years) (group 1A) and 45 healthy controls (11 males, 34 females; mean age 32.18±9.54; range, 21 to 57 years) (group 1B) to discover potential susceptibility variants. Subsequently, we verified these RA-susceptibility single nucleotide variants (SNVs) by genotyping an independent group of 500 Chinese subjects with RA (71 males, 429 females; mean age 48.6±12.2 years; range, 24 to 92 years) (group 2A) and 500 healthy controls (66 males, 434 females; mean age 32.3±10.1 years; range, 21 to 73 years) (group 2B). The study protocol was approved by the National Healthcare Group Ethics Committee (NHG DSRB 2014/01141). A written informed consent was obtained from each participant. The study was conducted in accordance with the principles of the Declaration of Helsinki.

We calculated the power of our study with 48 patients and 45 controls and significance level of 5%.[15] If the susceptibility variant has a prevalence of 5% in the controls and 20% in the patients, the power is 62.4%. If the variant has a prevalence of 5% in the controls and 30% in the patients, the power is then 92%. If the variant has a prevalence of 2% in the controls and 20% in the patients, the power is 83.1%.

Library preparation of the genomic deoxyribonucleic acid was performed with the NimbleGen SeqCap EZ kit (Roche, Penzberg, Germany). The targeted next-generation sequencing (NGS) panel was designed to capture the exons and the 5’ and 3’ untranslated region (UTR) up to 2kb upstream (Table 1). We sequenced the coding exons and flanking noncoding regions using Illumina MiSeq (San Diego, CA, USA). Amplification of the libraries was performed and assessed using a bioanalyzer (Agilent Technologies Inc., Santa Clara, CA, USA) and Qubit fluorometer (ThermoFisher Scientific, Waltham, MA, USA). The indexed samples were pooled at a final concentration of 8 pmol/L, and then parallel sequenced using MiSeq.

The resulting FASTQ files were aligned to the human reference genome (hg19), followed by duplicate removal by Picard Tools version 1.129. The Genome Analysis Toolkit version 3.3 was used for local realignment around indels, base recalibration, and variant recalibration and genotyping. Poor quality reads were removed and variants were selected based on quality score >30, coverage >10, and Hardy-Weinberg equilibrium (HWE) >0.001. Variants were annotated with dbSNP, UCSC, ClinVar, OMIM, KEGG, Pfam, Ensembl, ESP, 1000G, CADD, Polyphen, SIFT, ENCODE, HPRD, COSMIC, GERP, FitCons, and VISTA database using GEnome MINIing (GEMINI). The single-variant GWAS was performed using standard allelic testing implemented in PLINK tool. For multiple variant testing, the C-alpha variance-component test implemented in GEMINI was used with multiple allele frequencies thresholds.[16]

Candidate variants were prioritized through association analysis, selecting those with minor allele frequency (MAF) below 5% or predicted to be protein altering (i.e., missense, nonsense, frameshift or canonical splice site changes). Validation was performed in 500 cases and 500 controls using Sequenom’s MassARRAY system (Agena Biosciences, San Diego, CA, USA). One SNP, rs11567882 (GATA3), was genotyped using TaqMan (Applied Biosystems, Foster City, CA, USA) in the same set of 1,000 samples.

Statistical analysis

We calculated the power of our study with 48 patients and 45 controls and significance level of 5%.[15] If the susceptibility variant has a prevalence of 5% in the controls and 20% in the patients, the power is 62.4%. If the variant has a prevalence of 5% in the controls and 30% in the patients, the power is then 92%. If the variant has a prevalence of 2% in the controls and 20% in the patients, the power is 83.1%.

Results

The overall target enrichment and NGS yielded an average of 1,472,141 reads; 88% of these reads mapped to the targeted regions. The mean read length was 148.57 and GC content was 44.34%. After excluding variants that were off target, coverage lower than 10x and quality score below 30, we obtained a set of 3,512 variants. Among these, 2,193 (62.4%) were low-frequency variants with MAF <5%; of which 1,237 (35.2%) were novel to the 1,000 genome project dataset.

We used three strategies for identifying potential susceptibility variants. First, we examined 148 variants with significant association (p<0.05); 70 of them were low-frequency variants of MAF <5%. Fifty-three of the variants reside in the non-coding region, upstream of 5’UTR or 3’UTR. Two long non-coding ribonucleic acid variants and six synonymous variants were excluded. We further removed 32 HLA-DRB1 variants for subsequent analysis, except for rs9269688 (odds ratio [OR] 2.49), rs3180268 (OR 2.97), and rs77637983 (OR 0.43). Second, we conducted a pathway enrichment analysis using KEGG. Loss-of-function variants were identified in HLA-DRB1 rs9269957 and rs9269958 of the control group. We did not observe any other pathway for coding variants for RA in which a significant enrichment existed at p<0.05. Third, we investigated gene-based tests for the potential involvement of rare variants and RA susceptibility by using the C-alpha test. We identified an association between low-frequency variants and RA in five genes with p <0.05. There were nine non-synonymous low-frequency variants associated with RA. Other than PXK rs199881366, the other variants appeared to be protective. Thirty-nine low-frequency variants with significant association were identified for subsequent validation (Table 2).

In group 2A of 500 ACPA-positive RA patients, 85.5% were female, the baseline disease activity score-28 was 3.11±1.24, and 91.33% were positive for the rheumatoid factor. Group 2B consisted of healthy controls, of whom 13.2% were female.

Out of the 39 SNVs shortlisted for validation, 18 failed the MassARRAY multiplex assay design. Twenty-two SNVs were validated in 500 cases and controls, with 21 SNVs genotyped using iPLEX mass spectrometry (Agena Bioscience, San Diego, CA, USA) and one SNP using TaqMan chemistry (Table 3). On the MassARRAY system, results for three HLA-DRB1 SNVs and AHNAK2 were unsatisfactory due to cluster skewing and inefficient primer binding. Eight SNVs were found to be monomorphic. One SNV showed statistically significant association: the missense variant in IL-6 signal transducer (IL-6ST) rs777853685 (p.Cys47Phe) (p=0.0194).

Discussion

We sequenced 128 known RA susceptibility genes to determine if there are low-frequency variants within the coding and flanking noncoding regions and found one novel variant in IL-6ST. This gene was first identified as a susceptibility factor through meta-analysis.[17] Our variant is ~171kb away from the previously reported rs6859219.

The IL-6 is produced in the synovium of RA patients.[18,19] Variants in the IL-6 gene are associated with RA. IL-6ST is a 130 kDa signal-transducing component of the IL-6 receptor for IL-6. It is also involved in the signalling of ciliary neurotrophic factor, leukemia inhibitory factor, and oncostatin M. IL-6ST variant at position 5:55260065 (rs777853685) has a low population frequency of 0.3% in Genome Aggregation Database and multiple lines of computational evidence support a deleterious effect on gene product with Deleterious Annotation of Genetic Variants score of 0.8728. The alteration at this position represents a polar amino acid (cysteine) substitution with a non-polar hydrophobic residue (phenylalanine), which is likely to impact secondary protein structure. A pLI 0.998 score indicates that the gene is extremely intolerant to loss of function mutation. A gain-of-function variant of the other 80 kDa component of the IL-6 receptor is associated with RA.[20] How a loss-of-function mutation in IL-6ST contributes to the pathogenesis of RA is unknown.

The problem of the missing heritability of RA remains unresolved. Alternative explanations are that the rare variants lie outside the known susceptibility genes in the intronic regions, or that we have underestimated the contribution of the common variants (necessitating re-examination of the role of twin studies or re-calculation of the contribution of SNPs).[21] Beyond genetics, the pathogenesis of RA remains unresolved after decades of devoted research.[22]

A shortcoming of this study is the small number of patients in the discovery set. The strengths of this study are the single ethnicity of the research participants and the uniform ACPA status.

In conclusion, our results suggest that a mutation in IL-6ST 5:55260065 confers risk of RA in Chinese patients in Singapore.

The TTSH Rheumatoid Arthritis Study Group

Angela Li-Huan Chan, Grace Yin Lai Chan, Madelynn Tsu-Li Chan, Faith Li-Ann Chia, Hiok Hee Chng, Choon Guan Chua, Hwee Siew Howe, Ee Tzun Koh, Li Wearn Koh, Kok Ooi Kong, Weng Giap Law, Samuel Lee Shang Ming, Khai Pang Leong, Tsui Yee Lian, Xin Rong Lim, Jess Mung Ee Loh, Mona Manghani, Justina Wei Lynn Tan, Sze-Chin Tan, Teck Choon Tan, Claire Teo Min-Li, Bernard Yu-Hor Thong.

Table 1

List of genes containing rheumatoid arthritis susceptibility loci that are re-sequenced in this project

SN	Gene	CDS Length (base pairs)	SN	Gene	CDS Length (base pairs)
1	ABHD6	1014	65	IL6	639
2	ACOXL	1743	66	IL6R	1407
3	AFF3	3756	67	IL6ST	2757
4	AHNAK2	17388	68	ILF3	2697
5	AIRE	1638	69	INPP5B	2982
6	ANKRD55	1845	70	IRAK1	2139
7	ANXA3	972	71	IRF4	1356
8	ARAP1	4353	72	IRF5	1545
9	ARID5B	3567	73	IRF8	1281
10	ATG5	828	74	JAZF1	732
11	ATM	9171	75	LBH	336
12	B3GNT2	1194	76	LOC100506023	*
13	BLK	1518	77	LOC100506403	*
14	C1QBP	849	78	LOC145837	*
15	C5	5031	79	LOC339442	*
16	C5orf30	621	80	LY9	1968
17	CASP8	1617	81	MBP	915
18	CCL19	315	82	MED1	4746
19	CCL21	405	83	MMEL1	2340
20	CCR6	1125	84	MTF1	2262
21	CD2	1056	85	NFKBIE	1503
22	CD226	1011	86	P2RY10	1020
23	CD244	1113	87	PADI4	1992
24	CD28	705	88	PLCL2	3384
25	CD40	834	89	PLD4	1542
26	CD5	1488	90	PPIL4	1479
27	CD83	618	91	PRKCH	2052
28	CDK2	1041	92	PRKCQ	2121
29	CDK4	912	93	PTPN11	1782
30	CDK6	981	94	PTPN2	1248
31	CEP57	1503	95	PTPN22	2424
32	CFLAR	1443	96	PTPRC	3921
33	CLNK	1287	97	PVT1	*
34	COG6	1974	98	PYX	1737
35	CSF2	435	99	RAD51B	1278
36	CSF3	624	100	RAG1	3132
37	CTLA4	672	101	RAG2	1584
38	CXCR5	1119	102	RASGRP1	2547
39	DNASE1L3	918	103	RCAN1	834
40	EOMES	2118	104	REL	1860
41	ETS1	1458	105	RTKN2	1830
42	ETV7	1026	106	RUNX1	1443
43	FADS1	1335	107	SFTPD	1128
44	FADS2	1338	108	SH2B3	1728
45	FADS3	954	109	SMIM20	507
46	FCGR2A	933	110	SPRED2	1298
47	FCGR2B	810	111	STAT4	2247
48	FCGR3B	2229	112	SYNGR1	702
49	FCRL3	1335	113	TAGAP	2196
50	GATA3	1878	114	TEC	1896
51	GRHL2	810	115	TNF	702
52	HLA-DRB1	1422	116	TNFAIP3	2373
53	ICOSLG	1071	117	TNFREF14	852
54	IFNGR2	1530	118	TNFRSF9	768
55	IKZF3	537	119	TPD52	744
56	IL10	810	120	TRAF1	1251
57	IL1B	462	121	TRAF6	1569
58	IL2	936	122	TXNDC11	2958
59	IL20RB	489	123	TYK2	3564
60	IL21	489	124	UBASH3A	1986
61	IL23R	1890	125	UBE2L3	639
62	IL23RA	819	126	WDFY4	9555
63	IL2RB	1656	127	YDJC	972
64	IL3	459	128	ZNF438	2487
CDS: Coding sequence.

Table 2

Thirty-nine potential rheumatoid arthritis-associated variants detected in 21 genes through next-generation sequencing of 128 genes

Number	Gene	SNV	P	Chromosome	Reference base	Alternate base	Impact
1	C1orf141	1:67630314	0.02942	1	T	C	Intron variant
2	FCGR2A	rs1131184	0.0394	1	T	A	3 prime UTR variant
3	FCGR2A	rs12046367	0.002902	1	A	T	3 prime UTR variant
4	FCGR2B	1:161632071	0.04887	1	T	C	Intron variant
5	FCGR2B	rs375055702	0.03061	1	G	T	Intron variant
6	FCGR2B	rs56308545	0.03666	1	A	G	Intron variant
7	IL6R	rs76419864	0.0007269	1	T	C	Upstream gene variant
8	CD28	rs28541784	1.69E-07	2	C	T	Upstream gene variant
9	CD28	rs28688913	5.65E-07	2	C	T	Upstream gene variant
10	CD28	rs28718975	5.14E-06	2	T	C	Upstream gene variant
11	PXK*	3:58376372	-	3	C	A	Missense variant
12	PXK*	3:58382826	-	3	A	G	Missense variant
13	PXK*	rs199881366	-	3	T	C	Missense variant
14	EOMES	rs12715125	0.004245	3	C	G	Missense variant
15	TEC	4:48271895	0.04527	4	C	T	Upstream gene variant
16	TEC	4:48271899	0.03937	4	C	G	Upstream gene variant
17	IL6ST*	5:55237103	-	5	T	G	Missense variant
18	IL6ST*	5:55260065	-	5	C	A	Missense variant
19	IL6ST*	rs146973784	-	5	A	T	Missense variant
20	CCR6	6:167552230	0.02443	6	T	C	3 prime UTR variant
21	CCR6	6:167552236	0.0319	6	A	G	3 prime UTR variant
22	ETV7	rs369515633	0.04182	6	A	T	5 prime UTR variant
23	HLA-DRB1	rs3180268	0.001196	6	A	G	3 prime UTR variant
24	HLA-DRB1	rs77637983	0.01095	6	G	C	Missense variant
25	HLA-DRB1	rs9269688	0.006133	6	A	G	3 prime UTR variant
26	JAZF1	7:28220297	0.009422	7	G	A	5 prime UTR variant
27	IRF5*	7:128587373	-	7	C	T	Missense variant
28	IRF5*	rs113806178	-	7	G	A	Missense variant
29	IRF5*	rs201243166	-	7	C	T	Missense variant
30	C5	9:123814552	0.03006	9	G	T	Upstream gene variant
31	GATA3	rs11567882	0.04523	10	T	C	Upstream gene variant
32	GATA3	rs28395794	0.004603	10	A	G	5 prime UTR variant
33	CD5	11:60868006	0.0006823	11	A	C	Upstream gene variant
34	CD5	rs375347163	0.03061	11	C	G	3 prime UTR variant
35	CD5	rs72912997	0.04337	11	A	G	3 prime UTR variant
36	AHNAK2	rs2819428	0.04304	14	G	A	Missense variant
37	IRF8	rs75590645	0.023	16	G	A	3 prime UTR variant
38	IFNGR2	21:34775413	0.0432	21	G	A	5 prime UTR variant
39	UBE2L3	22:21920281	0.03778	22	T	A	Intron variant
SNV: Single nucleotide variant; UTR: Untranslated region; Asterisks indicate that variants were identified by C-alpha, while the rest was identified through association analysis.

Table 3

Genotype and allelic analysis of 10 of 22 single nucleotide variants in validation cohort

Assay	Gene	SNV	Genotype	Controls (n=500)	Cases (n=500)	Allelic frequency OR (95% CI)	p
Sequenom	ETV7	rs369515633	TT	341	341	1.04 (0.85-1.28)	0.6918
			TA	84	78
			AA	74	81
	FCGR2A	rs12046367	TT	362	340	1.22 (0.97-1.55)	0.0943
			TA	111	139
			AA	19	20
	GATA3	rs28395794	GG	436	449	0.84 (0.58-1.23)	0.3733
			GA	59	48
			AA	2	3
	IL6R	rs76419864	CC	500	498	2.51 (0.12-52.31)	0.5534
			CT	0	2
			TT	0	0
	IL6ST	5:55237103	TT	498	493	3.52 (0.73-16.97)	0.0948
			TG	2	7
			GG	0	0
	IL6ST	5:55260065	CC	499	492	8.06 (1.00-64.53)	0.0194
			CA	1	8
			AA	0	0
	IRF5	rs201243166	CC	496	497	0.75 (0.17-3.36)	0.7059
			CT	4	3
			TT	0	0
	IRF8	rs75590645	AA	452	445	1.198 (0.81-1.78)	0.3674
			AG	48	53
			GG	0	2
	PXK	rs199881366	TT	499	494	6.03 (0.72-50.18)	0.0965
			TC	1	6
			CC	0	0
TaqMan	GATA3	rs11567882	CC	439	446	0.865 (0.58-1.29)	0.4535
			CT	57	47
			TT	2	3
SNV: Single nucleotide variant; OR: Odds ratio; CI: Confidence interval; Of 22 single nucleotide variants chosen for validation, four were unsuccessful because of technical reasons and eight were monomorphic (data not shown).