The influence of polygenic risk scores on heritability of anti-CCP level in RA.

The objective of this study was to study genetic factors that influence quantitative anticyclic citrullinated peptide (anti-CCP) antibody levels in RA patients. We carried out a genome-wide association study (GWAS) meta-analysis using 1975 anti-CCP+ RA patients from three large cohorts, the Brigham Rheumatoid Arthritis Sequential Study (BRASS), North American Rheumatoid Arthritis Consortium (NARAC) and the Epidemiological Investigation of RA (EIRA). We also carried out a genome-wide complex trait analysis (GCTA) to estimate the heritability of anti-CCP levels. GWAS-meta-analysis showed that anti-CCP levels were most strongly associated with the human leukocyte antigen (HLA) region with a P-value of 2 × 10(-11) for rs1980493. There were 112 SNPs in this region that exceeded the genome-wide significance threshold of 5 × 10(-8), and all were in linkage disequilibrium (LD) with the HLA- DRB1*03 allele with LD r(2) in the range of 0.25-0.88. Suggestive novel associations outside of the HLA region were also observed for rs8063248 (near the GP2 gene) with a P-value of 3 × 10(-7). None of the known RA risk alleles (∼52 loci) were associated with anti-CCP level. Heritability analysis estimated that 44% of anti-CCP variation was attributable to genetic factors captured by GWAS variants. In summary, anti-CCP level is a heritable trait, and HLA-DR3 and GP2 are associated with lower anti-CCP levels.

INTRODUCTION

Anti-cyclic citrullinated peptide (anti-CCP) antibody is an important biomarker for rheumatoid arthritis (RA). It is more sensitive and specific than rheumatoid factor for diagnosing RA. Multiple studies have shown that presence of anti-CCP antibodies and level of anti-CCP in RA correlates with radiographic progression as well as extra-articular manifestations[1-5]. Higher anti-CCP levels are also associated with more active disease, more severe joint damage, worse functional disability and reduced quality of life[6, 7], and thus can be considered a proxy for disease severity. Anti-CCP antibody levels are relatively stable over time, have little correlation with disease duration, and change little after treatment[6], making anti-CCP a good target quantitative trait for genetic studies.

The anti-CCP test is a commercial assay that relies upon reactivity to several modified or synthetic cyclic citrullinated peptides. Research assays have been developed for autoantibodies directed against specific proteins and/or peptides found in the synovium such as vimentin, fibrinogen and enolase [8-13] in an effort to identify a specific antigen responsible for inciting RA. Multiplex assays for anti-citrullinated peptide antibodies (ACPA) demonstrate strong correlation between quantitative ACPA count and quantitative CCP levels, suggesting that the CCP assay is a valid proxy for more specific autoantibodies[12, 13].

RA heritability estimates suggest that genetic factors explain ~53-68% of RA susceptibility[14]. Since anti-CCP level is a measure of RA disease severity, genetic factors may also influence anti-CCP levels. Most previous studies of genetic determinants of anti-CCP have been limited to the human leukocyte antigen (HLA) region[15, 16] and a small number of RA candidate genes [17]. Most studies of the HLA-DRB1 Shared Epitope (SE), the strongest genetic risk factor for RA, have demonstrated consistent findings that the HLA-SE is associated with anti-CCP positivity; while some have found that the HLA-SE is associated with higher quantitative anti-CCP level specifically[15, 16]. In contrast to HLA-SE, Irigoyen et al reported that HLA- DRB1*03 was associated with lower anti-CCP levels in a Caucasian population[15], and Bang et al reported that HLA- DRB1*0901 was associated with lower anti-CCP levels in an Asian population[18]. However, the findings related to anti-CCP and other RA risk alleles are not convincing.

Genome wide association studies (GWAS) have been expanding for the past decade, and most RA genetic studies have focused on susceptibility analyses. To date, there are approximately 52+ loci validated as RA risk alleles outside the HLA-SE by GWAS and GWAS-meta analysis[19-23]. However, very few studies have focused on the genetics of RA disease severity, especially defined by anti-CCP level. Using Affymetrix 100K chip data from 531 patients in the Brigham Rheumatoid Arthritis Sequential Study (BRASS), we previously reported that HLA- DRB1*03 was associated with lower anti-CCP levels. We also found suggestive associations for SNPs in the HLA region that explained additional anti-CCP variation after adjusting for the DRB1*03 allele. These results were replicated in the North American Rheumatoid Arthritis Consortium (NARAC)[24].

In this study, we extend our GWAS on anti-CCP level to perform a meta-analysis including 3 large cohorts, BRASS, NARAC and the Epidemiological Investigation of RA (EIRA). Our purpose is to identify additional common genetic variants that may influence anti-CCP level to identify novel pathways associated with disease severity in RA. The current study provides unique information by extending the sample size to improve statistical power and using genotype data from denser and newer versions of the micro-array assay. In addition to GWAS analysis, we also apply two new methods, polygenic analysis and mixed linear modeling analysis, to assess the aggregated effect of common SNPS on anti-CCP level.

RESULTS

All subjects met 1987 ACR classification criteria for RA or were diagnosed by a board-certified rheumatologist. All subjects were anti-CCP positive, assessed by a second generation anti-CCP assay. 474 BRASS subjects, 823 NARAC subjects and 678 EIRA subjects passed all quality controls. The HLA region was the top ranked result with a p value of 2×10−11 (Figure 1) from meta-analysis. Table 1 shows the independent top 20 loci after removing the HLA region.

Figure 1

Genome-wide association analysis results for anti-CCP level.

Shown are strengths of association (−log10 P-value) for each SNP versus position along chromosomes 1 to 22.

Table 1

Top 20 loci associated with anti-CCP level outside the MHC region.

Chromosome	SNP	gene	bp	ma[a]	maf[b]	beta	p[c]	q[d]	Genotyped[e]
16	rs8063248	GP2	20166228	A	0.08	−0.29	2.9E-07	0.61	1
7	rs11766683	ELMO1	36801906	C	0.48	0.14	5.6E-06	0.44	0
2	rs155124	ITGA4	182016183	G	0.04	−0.42	1.0E-05	0.50	0
4	rs10489063	CD133	15926529	A	0.08	−0.25	1.2E-05	0.80	0
9	rs12345359	PIP5K1B	70777493	A	0.01	0.65	1.2E-05	0.92	0
17	rs8067785	ZNF18	11871706	G	0.15	−0.19	1.4E-05	0.46	1
2	rs4669784	NTSR2	11955419	C	0.38	0.15	1.8E-05	0.81	0
4	rs13143069		63108288	G	0.47	0.14	2.2E-05	0.54	0
3	rs4527407	SLITRK3	166693890	C	0.40	−0.14	2.5E-05	0.95	0
14	rs4903534	C14orf4	76551456	A	0.47	−0.13	2.8E-05	0.48	1
10	rs1147594	SLIT1	98955091	A	0.01	−0.57	3.0E-05	0.32	0
10	rs1907336		77824762	G	0.16	0.18	3.1E-05	0.89	1
2	rs1530875	B3GNT7	231907722	C	0.14	0.24	3.2E-05	0.41	0
5	rs2024141	CENTD3	140745184	A	0.08	0.25	3.2E-05	0.47	0
3	rs12485750	OSBPLIO	31392045	A	0.09	−0.24	3.4E-05	0.30	0
18	rs7227796		60653976	A	0.12	0.20	3.5E-05	0.80	0
1	rs662997	BCL9	145488076	A	0.03	0.39	3.7E-05	0.79	0
2	rs12615247		59040711	A	0.28	0.14	3.8E-05	0.58	1
11	rs1469170	ARHGAP20	110246224	A	0.30	−0.14	4.0E-05	0.88	1
9	rs12342611		109800925	C	0.48	0.13	4.3E-05	0.43	1

ma: minor allele

maf: minor allele frequency

p value is adjusted for 3 eigenvectors

Cochran’s Q tests for heterogeneity

1: genotyped in any of the three datasets

We found no evidence of systematic bias with λGC equal to 0.997. The Quantile-quantile (Q-Q) plot is shown in Figure 2a. The Q-Q plot removing the HLA region is shown in Figure 2b. A cluster of SNPs outside the HLA region on chromosome 16 were associated with anti-CCP level (Figure 3). Several SNPs, including the top SNP rs8063248, were genotyped in at least one cohort, and effect sizes were consistent among 3 cohorts (small Cochran’s Q-value, Supplementary Table 1). The top SNP rs8063248 was 63kb upstream from the GP2 gene (Glycoprotein 2 /Zymogen granule protein 2), which is an auto-antigen for inflammatory bowel disease[25] (IBD, Crohn’s disease), p=3×10−7.

Figure 2

Quantile-quantile (QQ) plots.

Shown are comparisons between observed −log10(p) and expected −log10(p). (a) QQ plot for all SNPs; (b) QQ plot after HLA region removed.

Figure 3

Regional association plot.

Showing strengths of association (−log10 P-value) versus position around rs8063248 along chromosome 16.

For the HLA region at (6p22.2 (26.03 Mb) to 6p21.32 (33.59Mb)), the regional association was plotted in Figure 4a. After controlling for the top SNP rs1980493, the association signal was no longer significant (Figure 4b), indicating that HLA associations with anti-CCP level are due to LD with the top SNP, or LD with causal variants in LD with the top SNP.

Figure 4

Association results for HLA region.

(a) Strengths of association (−log10 P-value) versus position along HLA region; (b) Association (−log10 P-value) versus position along HLA region adjusted for the top SNP rs1980493.

Using high resolution HLA typing in the NARAC samples, we explored LD between the individual HLA alleles, HLA-DRB1*0401, HLA-DRB1*0404 HLA-DRB1*03 etc., and top SNPs in the HLA region that were associated with anti-CCP level. Most of the top HLA SNPs were in LD with HLA- DRB1*03 with the highest r2 = 0.88 (Table 2). DRB1*03 was associated with low anti-CCP level in NARAC samples with p=0.0001, in BRASS samples with p=0.0003, in EIRA (imputed) samples with p=0.0003, with overall p=3×10−11. Table 3 shows the association between HLA-DRB1 subtype and anti-CCP level. Association between amino acid sites and anti-CCP also supports that HLA-DRB1*03 is associated with anti-CCP level (Supplementary Table 2).

Table 2

Linkage disequilibrium between top MHC region SNPs and HLA-DRB1 alleles.

rs	association with anti-CCP level p*	drb0101	drb0401	drb0404	drb03
rs2856674	2.39E-11	0.0001	0.02	0.006	0.88
rs9275576	2.86E-11	0.0001	0.02	0.006	0.86
rs2647044	3.19E-11	0.0001	0.02	0.006	0.86
rs1794520	4.20E-11	0.0001	0.01	0.006	0.85
rs3135353	3.58E-11	0.0004	0.02	0.008	0.71
rs3129861	2.69E-11	0.0009	0.02	0.008	0.69
rs2227138	3.66E-11	0.0004	0.02	0.008	0.69
rs1980493	2.23E-11	0.0009	0.02	0.01	0.64
rs3117097	3.23E-11	0.0009	0.02	0.01	0.64
rs3129953	4.71E-11	0.0009	0.02	0.01	0.64

p-value adjusted for 3 eigenvectors

Table 3

Association between HLA-DRB1 subtype and anti-CCP level.

HLA DRB1 type	frequency	beta	se	p*	q**
DRB1_03	0.08	−0.40	0.06	3.3E-11	0.59
DRB1_11	0.04	0.22	0.08	0.00976	0.31
DRB1_15	0.09	0.13	0.05	0.014	0.33
DRB1_01	0.15	0.10	0.05	0.026	0.11
DRB1_0101	0.13	0.13	0.05	0.010	0.13
DRB1_08	0.02	−0.24	0.12	0.052	0.89
DRB1_10	0.02	0.17	0.11	0.11	0.62
DRB1_07	0.05	0.09	0.07	0.20	0.63
DRB1_13	0.04	−0.10	0.08	0.24	0.01
DRB1_04	0.46	−0.02	0.03	0.56	0.42
DRB1_0401	0.31	0.04	0.04	0.26	0.57
DRB1_0404	0.10	−0.17	0.06	0.002	0.13
DRB1_0408	0.02	0.26	0.12	0.03	0.15
DRB1_09	0.01	−0.07	0.13	0.61	0.10

p value is adjusted for 3 eigenvectors

Cochran’s Q tests for heterogeneity

No established RA risk alleles[26] were associated with anti-CCP level from the meta-analysis after Bonferroni correction (Supplementary Table 3).

For the polygenic analysis, when we used EIRA as the discovery set, the polygenic risk score was significantly associated with anti-CCP level in BRASS. The highest variance explained was 1.7%, using a SNP set with a threshold of p<0.05. The association between the polygenic risk score and anti-CCP level was not significant in NARAC (Figure 5a). We found similar results using BRASS as a discovery set, namely that the polygenic risk score was associated with anti-CCP level in EIRA but not in NARAC (Figure 5b). Using NARAC as the discovery set, the polygenic risk score was not associated with anti-CCP level in either EIRA or BRASS (data not shown). GRS composed using RA risk alleles was not associated with anti-CCP level.

Figure 5

Association of polygenic risk scores with anti-CCP levels in the independent test datasets.

Variance explained by the PRS in the target samples for SNP set derived from discovery dataset for nine significance thresholds (p<10−5, 10−4, 10−3, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4). The y axis indicates R2. The significance level was indicated next to the threshold value. (a) EIRA as discovery set, BRASS, NARAC as test sets; (b) BRASS as discovery set, EIRA, NARAC as test sets. *: P-value ≤0.05; **: P-value ≤0.01.

In BRASS, EIRA and NARAC (the largest subset N=429 was used), genetic information summarized from the whole genome contributed to 65%, 37% and 41%, respectively, of total anti-CCP variation. The heritability estimates (Table 4) were not significant in BRASS, EIRA and NARAC (p=0.12, 0.15 and 0.2 respectively) due to the large standard error caused by small sample sizes. When we combined all 3 cohorts, heritability was estimated as 44% with a significant of p-value of 0.04.

Table 4

Heritability estimates for quantitative anti-CCP level.

Cohort	h2*	se	p
BRASS	0.65	0.54	0.12
EIRA	0.37	0.36	0.15
NARAC	0.41	0.58	0.2
BRASS, NARAC, EIRA	0.44	0.27	0.04

h2: heritability (0- 1)

DISCUSSION

We carried out a meta-analysis of GWAS for anti-CCP level from 3 cohorts comprising 1,975 subjects. We found that Anti-CCP level is a heritable trait. In addition to the HLA region, suggestive genetic effects outside the HLA region influenced anti-CCP level.

This study expands upon our previous 100K GWAS on anti-CCP titer for 531 BRASS subjects, followed by replication in NARAC only for the top 203 SNP findings. In the current report, we applied a different study design that used more homogeneous samples by restricting for anti-CCP+ cases, and excluding anti-CCP negative disease, which may have a different genetic architecture. We used a newer genotype platform with denser SNPs to achieve better power to identify genetic variants. We used imputation as a bridge to make BRASS, NARAC and EIRA samples comparable for the same set of SNPs. Since anti-CCP is not normally distributed, and measures from different assays are difficult to compare, we chose to use the inverse normal transformation method to rank the anti-CCP level. This normalization makes the estimated SNP effect comparable among cohorts, assuming that each cohort has similar anti-CCP level distributions, so that we can carry out a meta-analysis to combine GWAS results for all 3 cohorts. We also need to interpret the association with transformed anti-CCP level carefully. We would use the magnitude and the sign of beta to interpret the association, negative beta means the genetic variant is associated with lower CCP level, while the larger beta means a larger effect size.

Common HLA-SE risk alleles (DRB1*0101, DRB1*0401, DRB1*0404) were not strongly associated with anti-CCP level (p ~0.01 in NARAC, p ~0.002 for all 3 cohorts using imputed HLA type); however DRB1*03 was strongly associated with lower anti-CCP levels (p=0.0001 in NARAC, p=3×10−11 for all 3 cohorts), and explained 2.5% of the anti-CCP variance. The amino acid analysis provided additional support for the association between DRB1*03 (combination of 11S-13S-71K-74R[27]) and anti-CCP level. Specifically, the amino acid analysis identified HLA-DRB1 positions 11(Ser, S), 13(Ser, S), 71(Lys, K), 74(Arg, R) as the sites most significantly associated with anti-CCP level. The overall p value for these four amino acid sites are all significant (Supplementary table 2). This observation is consistent with a recent study of ~20,000 subjects, which reported that the DRB1*03 allele has a protective effect on RA, with a hazard ratio of 0.63[27] and with prior studies of DRB1*03 and anti-CCP level[15, 24]. LD between the top SNPs and the DRB1*03 allele are ~0.88 in NARAC and ~0.57 in BRASS. A regional plot of the results in the HLA region, before and after adjusting for the top SNP, suggests that association signals in the HLA region are due to one variant. However, due to the complexity of LD structure in the HLA region, it is challenging to identify the specific causal variant(s). Further strategies such as functional experiments are needed to discover the true causal gene/variant(s) in the HLA region.

In addition to the HLA region, the regional plot shows that a number of SNPs displayed suggestive associations with anti-CCP level on Chromosome 16 (Figure 3). Given the sample size, our study has limited power to detect moderate effect sizes to achieve genome wide significance. For example, our study had only 29% power to identify associations at the 5×10−8 significance level for the top chromosome 16 signal. We would have needed a sample size of 3131 to achieve 80% power. Also, GWAS assesses only common variants with minor allele frequency ≥ 5%. However, several criteria were applied to increase the likelihood of a true positive association. Specifically, we required that: 1) some of these SNPs were genotyped for at least one of the three cohorts, 2) the imputed SNPs had high accuracy scores (info score >0.8), and 3) association effect sizes were consistent among the three cohorts according to Cochran’s heterogeneity tests (Supplementary table 1).

The top finding outside HLA was on Chromosome 16, 3′ to the GP2 gene. GP2 plays an important role in bacterial and viral infections by recognizing antigen, thereby triggering antigen-specific immune responses[28]. In inflammatory bowel disease, antibodies to GP2 have been developed as diagnostic tools[25]. Since GP2 is immunogenic, it also may play a role in RA pathogenesis. Further studies of the role of GP2 in RA are necessary.

Literature searches focused on the remaining top SNPs revealed that ELMO1 is associated with diabetic nephropathy[29], whereas ITGA4 is a candidate gene for celiac disease[30]. CD133+ cells are endothelial progenitor cells that are mobilized in vasculogenesis in the inflammatory articular process that occurs in RA[31]. PIP5K1B is a candidate gene that may have a critical role in Sézary syndrome, a rare form of cutaneous T-cell lymphoma[32]. Although the relationship between these genes and anti-CCP is unknown, existing literature linking these genes and immunologic traits supports a role for these genes in RA.

Most prior genetic studies for anti-CCP were focused on anti-CCP status (+/−), showing an association between some RA risk alleles and anti-CCP status[33, 34]. Besides the reported association of HLA- DRB1*03 with lower anti-CCP level[15], one paper studied genetic associations with anti-CCP level and reported that one PADI4 haplotype, carrying a PADI4 risk allele, was associated with higher anti-CCP level in new-onset RA[17]. We examined 52 recently validated RA risk alleles[26] in the present study and found no evidence for association with anti-CCP level, suggesting that anti-CCP genetic architecture may differ from disease susceptibility architecture.

Our study does have some limitations. First, our analysis was limited to a study of the CCP test, a commercial assay that tests reactivity to several modified or synthetic cyclic citrullinated peptides[35], as opposed to ACPA assays, which test reactivity to individual peptides (e.g., vimentin, fibrinogen, enolase). Serum anti-CCP antibody levels are, however, highly correlated with total serum ACPA responses[9, 13, 36, 37] .

Second, although anti-CCP levels were measured in a single laboratory for BRASS and EIRA, anti-CCP was measured in three laboratories for NARAC. To examine the potential effects of differing techniques between these laboratories, we performed sensitivity analyses treating each NARAC laboratory subset as a separate cohort and transformed anti-CCP level separately for each subset. GWAS results from the meta-analysis combining data from BRASS, EIRA and the three NARAC sub-datasets were similar to results from the primary analysis, which utilized combined NARAC data (Supplemental figure 1). Because of small sample sizes, we did not use the NARAC subsets as discovery sets for the polygenic model. Because the heritability estimate is sensitive to the homogeneity of phenotype, heritability was estimated only for the largest NARAC sub-datasets (N=429), not the smaller two sub-datasets (N=154, 245 respectively).

In summary, we estimated that heritability for anti-CCP is approximately 44%. Given the modest sample size, the heritability estimates was not highly significant. However, the estimate was similar across cohorts, supporting the veracity of our findings. Independent samples will be needed to confirm our findings. Using the GCTA approach, we estimated that all SNPs from the MHC region explain ~3% of the anti-CCP variance, whereas the polygenic risk score explains an additional ~2% of the variance. Thus, much of the heritability remains unexplained. The challenge is to identify rare variants with larger effect sizes through next generation sequencing and to find additional common variants with low to moderate effect size using cohorts with larger sample sizes.

Conclusion

We confirmed that HLA is associated with anti-CCP level in our quantitative trait analysis. Specifically, HLA-DR3 and GP2 are associated with lower anti-CCP levels. Further studies of the functional role of HLA-DR3 and GP2 may provide insights into the mechanisms of disease severity and ultimately lead to targeted therapy if SNPs can predict patients who need more or less intense therapy. Although SNPs outside the HLA region also show modest associations with anti-CCP level, these associations do not reach genome wide levels of significance. Anti-CCP level is a heritable trait in which genetic factors explain ~44% variance based on heritability analysis, where HLA-DR3 explains ~2.5% of the variance. Based on polygenic results, the additive genetic effect of common alleles outside the HLA region only explains ~2% of the variance. Therefore, much of the heritability is ‘hidden’ and other techniques such as next generation sequencing should be pursued.

METHODS

Subjects

is a prospective, observational, hospital-based RA cohort established to identify biomarkers and genetic predictors of drug response and disease activity[38]. Patients were eligible if they had a diagnosis of RA and were older than 18 years. All study protocols were approved by the Brigham and Women’s Hospital Institutional Review Board and informed consent was obtained from all subjects. The cohort is predominantly female (82.3%) with a mean age of 58.0 (±13.6) years and average disease duration of 16.1 (±12.5) years.

samples were collected by 4 different study designs. “Family collection” samples were from multiplex families (primarily affected sibling pairs) in which at least one sibling had documented joint erosions. The other collections included samples from three prospective observational cohorts, the National Data Bank of Rheumatic Diseases, the National Inception Cohort of Rheumatoid Arthritis, and the Study of New Onset Rheumatoid Arthritis. A more detailed description can be found elsewhere[33]. 74.1% of the cohort were female, with a mean age of 56.0 (±12.3) years and average disease duration of 10.6 (11.1) years.

is a population based case control study carried out in Sweden, comprising residents of south and central Sweden between the ages of 18 and 70 years in 1996 to 2005[39]. An incident case is defined as a person in the collection who for the first time receives a diagnosis of RA (newly diagnosed cases) according to the 1987 ACR criteria. 72.3% of EIRA cases were females, with a mean age of 50.8(±12.8) years, and all are new onset cases (<2 years from the onset of sympotoms, with 85% within 1 year after the first onset of symptoms).

Outcome

Anti-CCP levels were measured using second generation ELISA assays (NARAC and BRASS: Inova Diagnostics; EIRA: Immunoscan-RA Mark2 ELISA). Anti-CCP antibody levels were transformed to a normal distribution by taking the inverse normal of the rank in BRASS, NARAC and EIRA separately. (Supplementary Figure 2 shows the anti-CCP distribution before and after transformation).

Genotyping

BRASS samples were genotyped using the Affymetrix 6.0 chip (900K). NARAC samples were genotyped using the Illumina 550K chip. EIRA samples were genotyped using the Illumina 317K chip. A quality control (QC) process was performed in these 3 datasets separately, where QC steps included filtering SNPs with >1% missing data and individuals with >5% missing data, followed by filtering SNPs with minor allele frequency (MAF) <1% and a Chi-Square test of Hardy Weinberg equilibrium PHWE<10−5. We then used individual pair-wise identity-by-state estimates to remove related and potentially contaminated samples. Data processing and QC were performed in PLINK [40]. To address population stratification and remove outliers, Principle Components Analysis (PCA) was done on the combined dataset, using more than 30K high quality, commonly genotyped SNPs across three cohorts, and genetic outliers were removed using the default setting in EIGENSTRAT[41] (Supplementary Figure 3 shows the PCA result after the outliers were removed).

Imputation

Imputation was conducted on genotype data for each data set separately from IMPUTE software, using haplotype-phased HapMap Phase 2 European CEU founders as the reference panel, allowed us to merge data across different genotype platforms[42]. Imputation yielded posterior genotype probabilities as well as imputation quality scores at SNPs not genotyped with a minor allele frequency ≥1% in the HapMap CEU. We removed SNPs with imputation quality scores <0.5 from our final analysis. Imputation for HLA type was also conducted in SNP2HLA software[43], using the Type 1 Diabetes Genetics Consortium collection as high quality HLA reference panel (N=>5,000)[43].

GWAS

Using additive models, association analyses were performed for more than 2 million genotyped/imputed SNPs with transformed anti-CCP values in each dataset. To adjust for residual population stratification effects, the first three eigenvectors from the PCA were controlled for in each model. All association analyses were done using SNPTEST[42]. Inverse variance weighted meta-analysis was utilized to combine GWAS results for all 3 cohorts to identify loci associated with anti-CCP level (N=1,975). We also conducted Cochran’s Q tests for heterogeneity across collections using the β coefficients. We did not consider SNPs if they showed strong evidence of heterogeneity across the collections (Cochran’s Q P<0.001). Quantile-quantile (Q-Q) plots and genomic control inflation factor-λGC were employed to illustrate any unmeasured systematic bias. We studied 52 non-HLA risk alleles for RA for association with CCP level based on a recently published meta-analysis[26].

Polygenic modeling

Since genome-wide significance criteria are stringent, most genetic variations will not achieve those thresholds. Polygenic modeling is a method to evaluate the aggregate effects of many variants with small effects. Information from different SNPs are combined into a ‘risk’ score to test whether that score has any association with anti-CCP level[44]. BRASS, NARAC and EIRA were each used as discovery sets using the same GWAS approach, adjusted for the PCA. Linkage disequilibrium (LD) pruning (LD r2>0.1) was done with respect to the association p-value to obtain a set of independent SNPs, which were maximally associated with anti-CCP level. The PRS then was computed in each of the other two test sets. SNP sets were selected for nine different threshold p-values (p<10−5, 10−4, 10−3, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4). For each SNP set, beta-weighted allele counts were summed into a polygenic risk score. The polygenic risk score was calculated for each validation set individual i as where is the discovery GWAS β- coefficient with respect to anti-CCP level, and dij is the individual i’s posterior probability dosage (from zero to two) of that allele. The resulting PRS was tested for association with anti-CCP level using linear regression modeling, adjusting for 3 eigenvectors. P-values for the association between polygenic risk score and anti-CCP level, the numbers of SNPs and estimated variance explained were reported. Specific PRS (genetic risk score, GRS) was also composed using RA risk alleles. Association between GRS and anti-CCP level was conducted.

Mixed linear model

Heritability was directly estimated from the BRASS, NARAC and EIRA GWAS datasets. Genome genotype data and imputed SNP data were summarized into a Kinship matrix described by Yang et al [45], where Aij is the relationship estimate for subject i and j for i≠j. Genetic variance then was estimated as the random effect coefficient of the kinship matrix in a MIXED model using the genome-wide complex trait analysis (GCTA) software[45]. The heritability estimate variance was calculated from the inverse of the variance-covariance matrix, and significance was assessed by the Wald test.