Genetic Association and Expression Correlation between Colony-Stimulating Factor 1 Gene Encoding M-CSF and Adult-Onset Still's Disease

Adult-onset Still's disease (AOSD) is a rare and inflammatory disorder characterized by spiking fever, rash, arthritis, and multisystemic involvement. HLA has been shown to be associated with AOSD; however, it could not explain the innate immunity and autoinflammatory characteristics of AOSD. To assess the genetic susceptibility of AOSD, we conducted a genome-wide association study (GWAS) on a cohort of 70 AOSD cases and 688 controls following a replication study of 36 cases and 200 controls and meta-analysis. The plasma concentrations of associated gene product were determined. The GWAS, replication, and combined sample analysis confirmed that SNP rs11102024 on 5′-upstream of CSF1 encoding macrophage colony-stimulating factor (M-CSF) was associated with AOSD (P = 1.20 × 10−8, OR (95% CI): 3.28 (2.25~4.79)). Plasma levels of M-CSF increased in AOSD patients (n = 82, median: 9.31 pg/mL), particularly in the cases with activity score ≥ 6 (n = 42, 10.94 pg/mL), compared to the healthy donors (n = 68, 5.31 pg/mL) (P < 0.0001). Patients carrying rs11102024TT genotype had higher M-CSF levels (median: 20.28 pg/mL) than those with AA genotype (6.82 pg/mL) (P < 0.0001) or AT genotype (11.61 pg/mL) (P = 0.027). Patients with systemic pattern outcome were associated with elevated M-CSF and frequently observed in TT carriers. Our data suggest that genetic variants near CSF1 are associated with AOSD and the rs11102024 T allele links to higher M-CSF levels and systemic outcome. These results provide a promising initiative for the early intervention and therapeutic target of AOSD. Further investigation is needed to have better understandings and the clinical implementation of genetic variants nearby CSF1 in AOSD.


Introduction
Adult-onset Still's disease (AOSD) is an inflammatory disorder characterized by spiking fever, macular rash, leukocytosis, arthritis, variable multisystemic involvement, and increase of acute phase reactants [1,2]. It is a rare disease with crude prevalence of only 1-2 cases in 100,000-1 million annually [3,4]. The clinical features and disease progression of AOSD vary considerably. In severe cases, AOSD may lend to permanent joint destruction, organomegaly, lymphadenopathy, serositis, and aseptic meningitis [5][6][7]. Due to its characteristics and predominate dysregulation of innate immunity, AOSD has been considered an autoinflammatory disease [8,9]. Standard treatments for AOSD include corticosteroids as the first-line, the use of nonsteroidal anti-inflammatory drugs (NSAIDs), and disease-modifying antirheumatic drugs (DMARDS) to manage the clinical symptoms [6,10,11]. Nevertheless, AOSD still lacks effective therapeutics, as its etiology and pathophysiology remain unclear [12].
Most of the previous genetic studies on AOSD revealed association with variants on human leukocyte antigen (HLA) class I and II regions, such as HLA-BW35, HLA-DRB1 * 15 and HLA-DRB1 * 04 [19,20]. However, the association was inconsistent and controversial among the various studies of different populations [27,28]. In addition, the results of genetic studies did not link to the pathogenesis, and none of them have been associated with AOSD disease outcomes [14,28]. Here, we enrolled 106 patients and 888 population controls and applied GWAS discovery and subsequent replication analysis to investigate the genetic susceptibility of AOSD. We explored the correlation between the identified genetic risk factor(s) and disease severity or outcomes and examined the functional implication of the associated genetic variants in AOSD.

Patient and Public Involvement.
This study was carried out following the rules of the Declaration of Helsinki of 1975, which was revised in 2013. This study was approved by the ethics committee of the Institutional Review Board (IRB) of Taichung Veterans General Hospital (CF13321), and the written consent was obtained from each participant. We enrolled a total of 106 AOSD patients fulfilling the Yamaguchi criteria [29] between January 2010 and December 2015. Patients with infections, malignancies, or other rheumatic diseases were excluded. The disease activity of AOSD was assessed with a modified Pouchot score described by Rau et al. [23]. According to the proposed classification of disease courses of AOSD [9], the AOSD patients with follow-up at least one year were classified into two patterns of disease outcomes: (i) the "systemic pattern" that includes the monocyclic and the polycyclic form and (ii) the "chronic articular pattern" (persistent arthritis involving at least one joint destruction and lasting longer than 6 months) [30]. All of the AOSD patients were unrelated Han Chinese. Seventy of them were randomly selected as the case group in the GWAS discovery cohort.
We recruited 924 ethnically and geographically matched healthy subjects as the population controls from a biobank under a nationwide population study, which comprises 9,980 Han Chinese descendants [31]. There was no selfreport of rheumatic diseases among the recruited controls from Taiwan, where 98% of the population is made up of Han Chinese. Of the 924 population controls, we randomly obtained 724 controls for the GWAS discovery cohort and the rest of 200 individuals were for the replication cohort, which was an independent analysis conducted by 36 AOSD cases versus 200 population controls to validate the statistically significant SNPs derived from the initial GWAS results.

Genotyping and Quality Controls in the Genome-Wide
Scan. Genomic DNA was extracted from the peripheral blood of the enrolled subjects using Flexi Gene DNA kits (Qiagen, Hilden, Germany). We performed GWAS on the samples obtained from 70 AOSD cases and 724 population controls using the Affymetrix SNPs Array 6.0 platform (Santa Clara, CA, USA), which is composed of 909,622 SNPs [31]. Briefly, 200 ng of genomic DNA of each sample was PCR amplified, fragmented, precipitated, and resuspended in the appropriate hybridization buffer. After hybridization, the BeadChip oligonucleotides were extended by a singlelabeled base, which were detected by fluorescence imaging with an Affymetrix Bead Array Reader. Normalized bead intensity data obtained for each sample were loaded into the Affymetrix SNP Array 6.0 software, which converts the fluorescence intensities into SNP genotypes. The genotype calls were generated using the Birdseed method (Birdseed v2) with Affymetrix Power Tools (version, apt-1.10.2).
We analyzed the GWAS data by the software Plink (v1.90b5) using logistic regression modeling with covariates: sex, and ancestry-specific principal components (i.e., PC1, PC2, PC3, PC4, PC5, PC6, PC7, PC8, PC9, and PC10). The genomic coordinates are based on NCBI Human Genome Build 37 (GRCh37). Of the 909,622 genotyped SNPs, 868,494 SNPs located on the autosomal chromosome were included for the quality control (QC) process. Our quality control for SNPs was a call rate higher than 0.95, and 59,494 variants were removed due to missing genotype data. We excluded 163,017 variants due to minor allele threshold(s) (minor allele frequency ðMAFÞ < 0:01). The quantile-quantile plot was generated using the R/Bioconductor package: GWASTools [32]. The genomic inflation factor was calculated by Plink (v1.90b5). Linkage disequilibrium (LD) analysis of the SNP rs11102024 was performed by using the LDproxy module of the online software package LDLink (https://analysistools.nci.nih.gov/LDlink).

Replication Analysis and Targeted Gene Sequencing.
For the replication study, we applied TaqMan assays (Thermo Fisher Scientific, CA, USA) or direct sequencing on the associated SNPs revealed by the initial GWAS using an independent sample set of 36 AOSD cases and 200 population controls. The oligonucleotide primers for polymerase chain reaction (PCR) amplification of the SNP rs11102024 near CSF1 gene are forward primer: 5 ′ -TCCTATTGCATTGG GCATATT-3 ′ , and reverse primer: 5 ′ -TCCATTTACGCC TCAACTCA-3 ′ , and those for SNP rs9636107 near TCF4 gene are forward primer: 5 ′ -GCTGGTACCAAGGAAAG CTG-3 ′ , and reverse primer: 5 ′ -CCTGCTGGTGTTTTGT TTTG-3′. The PCR reaction was performed in three steps: 3 min at 95°C; then 40 cycles of 20 seconds at 95°C, 30 seconds at 58°C, and 30 seconds at 72°C, followed by 7 minutes at 72°C.

Determination of Plasma Levels of Macrophage Colony-
Stimulating Factor (M-CSF). The plasma samples of 82 patients with AOSD were collected at the active status of the disease, and the patients received corticosteroids and/or the nonsteroidal anti-inflammatory drugs (NSAIDs). Besides, the disease-modifying antirheumatic drugs (DMARDs) had also been prescribed for the patients, which included methotrexate (patients number ðnÞ = 74), hydroxychloroquine (n = 66), cyclosporine (n = 30), sulfasalazine (n = 18), and azathioprine (n = 10). Levels of plasma M-CSF were determined on the samples from these 82 active AOSD patients and 68 population controls using enzyme-linked immunosorbent assay (ELISA) according to the manufacturer's instructions (Ray-Biotech Inc., GA, USA).

Statistical Analysis.
For the GWAS, we conducted the statistical analysis for the associations by comparing the allele frequencies between AOSD cases and population controls. The SNP association was examined by Fisher's exact tests and rank-ordered according to the lowest P value. All the P values were two-tailed. The corrected Pc values were adjusted using Bonferroni's correction for multiple comparisons (645,983 for GWAS SNPs), and Pc value for genomewide significance should be <7:8 × 10 -8 (0.05/645,983). The odds ratios (OR) were calculated using Haldane's modification [33]. Levels of plasma M-CSF in the different groups were compared by the nonparametric Mann-Whitney U test. The differences in the frequencies of significant alleles among AOSD patients with different courses of disease outcomes were examined using Fisher's exact test. The correlation coefficient was obtained by nonparametric Spearman's rank correlation test. A probability of less than 0.05 was considered to be significant.
3 Journal of Immunology Research population controls) passed the filters. The mean call rate is 98.3%. The SNPs on sex chromosomes were excluded from the GWAS analysis. The principal component analysis (PCA) map showed that there is no difference in the distribution of ancestry among AOSD cases and population controls (Supplemental Figure 1). A quantile-quantile plot of the test statistics was used for quality control and revealed that the population matching was successful (Supplemental Figure 2). The genomic inflation factor (λGC) [34] is 1.05, suggesting that the population structure of our GWAS is generally acceptable.
When comparing the allele frequencies of the 645,983 SNPs of 70 AOSD cases and 688 population controls, no SNP reached the threshold of genome-wide significance (Pc < 7:8 × 10 -8 ). Given AOSD is a rare disease and small sample size used in this study, we set the cut off P value as 1 × 10 -6 and the top 4 significant SNPs with highly suggestive association are shown in Figure 1 and Table 2. Of these 4 SNPs, two SNPs (rs35910146 and rs6948305) are on the upstream of pseudogenes, and rs9636107 does not match Hardy-Weinberg equilibrium (P = 0:03) in the dataset of 688 population controls ( Table 2). The SNP rs11102024 on 5 ′ -upstream of CSF1 (colony-stimulating factor 1) gene encoding for M-CSF (macrophage colony-stimulating factor) on chromosome 1p13 showed significant association with AOSD (P = 3:70 × 10 -7 , OR (95% CI): 3.27 (2.07~5.17)) ( Table 2). By comparison, the association strength between SNPs on HLA regions and AOSD was weaker than that of variants near CSF1 gene (Supplemental Table 1).
3.3. Replication, Meta-Analysis, and Linkage Disequilibrium of Variants near CSF1 Gene. We used the samples of an independent cohort (36 AOSD cases and 200 controls) to replicate the association. Among the 4 SNPs discovered by GWAS, only SNP rs11102024 near CSF1 displayed signifi-cant association with AOSD (P = 0:022, OR ð95%CIÞ = 2:47 (1.20~5.08)) ( Table 3). The meta-analysis of the two datasets from the initial GWAS and replication revealed that the P value of the heterogeneity test between studies (P het ) was 0.36, suggesting that there was no difference between both studies (P value of the heterogeneity test > 0:05) (Figure 2, Table 3). Subsequently, we performed combined sample analysis using the genotyping data of 106 AOSD cases and 888 population controls, and the SNP rs11102024 showed strong association with AOSD (P = 1:20 × 10 -8 ; OR (95% CI): 3.28 (2.25~4.79)) (Table 3, Figure 2). We used the genetic dataset of East Asian population of the 1,000 genome project for linkage disequilibrium (LD) analysis and identified 11 SNPs having strong LD with rs11102024 (r 2 > 0:7) (Supplemental Table 2, Supplemental Figure 3). We also sequenced the CSF1 gene; however, there were no missense SNPs with significant association with AOSD (data not shown).

Increased Plasma M-CSF Levels in AOSD Patients and
Correlation with the Activity Score of the Disease. We hypothesized that genetic variant near CSF1 might have an impact on the disease progress/outcome of AOSD. We determined the plasma levels of M-CSF in 82 patients who had available plasma samples. Significantly higher levels of plasma M-CSF were detected in AOSD patients (median 9.31 pg/mL, inter-quartile range (IQR) 6.21-16.91 pg/mL) compared with the healthy controls (n = 68) (median 5.31 pg/mL, IQR 4.12-6.85 pg/mL) (nonparametric Mann-Whitney U test, P < 0:0001) (Figure 3(a)). Moreover, significantly higher M-CSF levels were detected in the AOSD patients with activity score ≥ 6 (n = 42; median 10.94 pg/mL; IQR 6.73-20.23 pg/mL), compared to those with activity score between 3 and 5 (n = 40; 7.68 pg/mL, IQR 5.88-11.38 pg/mL) (Figure 3(b)).
Regarding the disease outcome, 64 of 82 AOSD patients had the systemic pattern, and 18 cases had the chronic articular pattern. Significantly higher M-CSF levels were observed in AOSD patients with the systemic pattern (median 10.79 pg/mL, IQR 6.40-18.53 pg/mL) compared to those with the chronic articular pattern (6.23 pg/mL, IQR 4.37-8.41 pg/mL) (P = 0:0007) (Figure 3(e)). In addition, the

Discussion
AOSD possesses similar clinical presentations with rheumatoid arthritis (RA); however, AOSD does not have RF or significant memorial lymphocyte involvement, and the genetic susceptibility revealed by this study showing no similarity with that of rheumatoid arthritis [35]. Consistent with the previous study, we did not find association between AOSD and MEFV gene mutations [36][37][38] or other genetic mutations of monogenic autoinflammatory disorders [39,40]. The systemic juvenile idiopathic arthritis (juvenile-onset Still's disease) has multiple-gene predisposition [8,41], including HLA-DRB1 [21,[42][43][44]. In addition, both HLA class I and II regions have been reported to be risk loci for AOSD [29]. However, these findings are controversial to the clinical characteristics of juvenile-onset or adult-onset Still's disease, as HLA would implicate the predominant adaptive immune response. In this study, we find that SNPs on HLA class I and class II had only weak association with AOSD, but 5′-upstream variants of CSF1 reached genomewide significance. The different observations may relate to the enrolling criteria for patients with diverse clinical presentations of Still's disease [18,44]. Further studies are needed to verify these genetic associations with AOSD in different populations. CSF1 gene encoding M-CSF has been implicated in the differentiation and activation of monocyte/macrophage lineage [25,26]. The increased M-CSF plasma levels have been considered a biomarker for AOSD [22,25,26]. In our study, M-CSF expression correlating with AOSD activity and severity is consistent with the results of previous studies [22,25,26]. Particularly, our GWAS first demonstrates that SNP rs11102024 on 5′-upstream of CSF1 reached genomewide significance, and the subsequently replication study confirmed the strong genetic association with AOSD. We further identified 11 SNPs having strong LD with rs11102024 (r 2 > 0: 7), and most of them show regulatory potential [45]. As the associated SNP rs11102024 is only 22 kb upstream of CSF1 and the nearby LD variants show regulatory potential, we propose that these SNPs may upregulate the expression levels of M-CSF, resulting in the disease of AOSD. In support of the functional roles of associated genetic variants, we indeed found significantly higher levels of plasma M-CSF in our AOSD patients with rs11102024 TT genotype. Taken together, our data suggest that the variants near CSF1 may play important pathogenesis roles to regulate the gene expression and participate in the disease course of AOSD. Further studies are needed to investigate the potential mechanism of action of the genetic variants in the context of AOSD.
The clinical course of AOSD can be divided into two main patterns with different prognoses: systemic pattern (monophasic or polycyclic) and chronic articular pattern. The systemic pattern is characterized by predominantly systemic features including fever, rash, serositis, and organomegaly with or without articular symptomatology. By comparison, the chronic articular pattern is characterized by the severe articular manifestations mimicking rheumatoid arthritis. In this study, we observed higher plasma levels of M-CSF in AOSD patients whose activity score ≥ 6, cases with systemic pattern, or carriers with rs11102024 TT genotype. These data suggest that rs11102024 T allele, which is highly associated with M-CSF plasma levels, could be used to predict the severity and system symptoms in the early phase of AOSD and act as a promising genetic marker for early intervention to improve AOSD outcome.
As AOSD is a polygenetic autoinflammatory disorder, increasing biologic agents are investigated to target its proinflammatory cytokines, such as IL-1 family (particularly IL-1β and IL-18), IL-6, and TNF-α. Growing evidences and clinical trials indicate that anticytokine biologic agents, e.g., anakinra (interleukin-1 receptor antagonist) [46] and tocilizumab (anti-IL-6 receptor antibody) [47], are becoming plausible therapeutic options for the management of AOSD. Recently, Jaguin et al. reported the phenotypic and genomic markers of M-CSF-generated human macrophages polarized towards M1 or M2 subtype upon the action of lipopolysaccharide and interferon-γ (for M1) or interleukin-(IL-) 4 (for M2) [48]. The ability of human M-CSF-generated macrophages to polarize toward M1 or M2 subtype was associated with enhanced secretion of TNF-α, IL-1β, IL-12p40, CXCL10, and IL-10 (for M1) or CCL22 (for M2) [48]. In addition, Bellora et al. reported that M-CSF could induce the expression of membrane-bound IL-18 in human blood monocytes differentiating toward macophages [49]. The enhancement of these macrophage-triggered cytokines was frequently observed in AOSD patients [22,50]. This functional property of M-CSF highlights the innate immunity and subsequent adaptive immune responses implicated in AOSD pathogenesis. In this study, we report that the novel SNP rs11102024   Journal of Immunology Research near CSF1 gene encoding for M-CSF is highly associated with AOSD and links to higher plasma levels of M-CSF, which would be a potential and promising therapeutic target for AOSD.

Conclusion
Adult-onset Still's disease (AOSD) is a rare, but systemic inflammatory disorder characterized by spiking fever, rash, These results are expressed as the mean ± SEM with each dot representing the data of an individual. Statistical analysis was performed using a two-tailed, nonparametric Mann-Whitney U test. 8 Journal of Immunology Research leukocytosis, arthritis, and multisystemic involvement, with increased macrophages and M-CSF. Herein, we identify a novel SNP rs11102024 on the 5 ′ -upstream of CSF1 gene to be significantly associated with AOSD in a GWAS and subsequent replication analysis. Reflecting the function of CSF1, encoding M-CSF which is involved in macrophage differentiation and inflammatory responses, the SNP rs11102024 demonstrated the significant association with plasma M-CSF levels, disease activity, and disease outcome of AOSD. These data suggest that the gene variant rs11102024 nearby CSF1 could be a potential prognostic factor for a disease outcome in AOSD. These results provide a promising initiative to predict and early intervention to improve the treatment and healthcare for AOSD. Further investigation is needed to a better understanding of the pathogenesis of AOSD and clinical implementation for the genetic variants of CSF1.

Data Availability
All data relevant to the study are included in the article or uploaded as supplementary information.

Conflicts of Interest
The authors declare no conflict of interest.

Authors' Contributions
All authors made substantive intellectual contributions to the present study and approved the final manuscript. Y-MC and W-TH contributed equally to this work, designed the study, conceived the study, performed clinical assessment as well as data acquisition, and drafted the manuscript. C-WH performed clinical assessments on study subjects, conducted the analysis of data, and drafted the manuscript. W-HC and J-LL designed the study and performed the data analysis. N-RG, W-CC, and Y-SL conceived the study, obtained the laboratory and genomic data, and performed the analysis. D-YC and S-I H contributed equally to this work, generated the original hypothesis, designed the study, acquired the clinical data and data analysis, and drafted and revised the manuscript.