Method of identifying a gene associated with a disease or pathological condition of the disease

Sample Collections

1. Patients

58 patients diagnosed as having RA were unrelated individuals of Han Chinese descent recruited from hospitals in Southern and Eastern China (Guangzhou and Changzhou) using 2010 Rheumatoid Arthritis Classification Criteria established by American College of Rheumatology and European League Against Rheumatism Collaborative Initiative (2010 ACR/EULAR).

In addition, 66 healthy and unrelated blood donors of Han Chinese ancestry from Medical Center for Physical Examination and Health Assessment, were included as controls.

Detailed descriptions of sequenced individuals and clinical characteristics of the enrolled patients are provided in Table 1 and 2. Written informed consent was obtained from all of the participants, and the study was registered in Chinese Clinical Trial Registry (ChiCTR-ROC-17010351) and approved by the local ethics committees of Macau University of Science and Technology (Macau, China).

2. Confirmation of the Ethnicity of the RA and Control Groups

The ethnicity of the patients suffering from RA, i.e. RA group, and the healthy control group was verified by conducting ancestry composition analysis using admixture v1.3.0 (https://www.genetics.ucla.edu/software/admixture) and multidimensional scaling in PLINK v1.07 (http://zzz.bwh.harvard.edu//plink/). The results are shown in FIG. 1 . Three ethnic populations were used as reference samples from 1000 Genome Project Phase III data (http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chr9.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz), including Utah Residents with Northern and Western European Ancestry (EUR-CEU), Yoruba in Ibadan, Nigeria (AFR-YRI) and Southern Han Chinese (EAS-CHS).

3. Preparation of the Samples for Gene Sequencing

Blood samples were collected from the patients of the RA group and the healthy people from the control group, according to protocols approved by local institutional review boards. Genomic DNA was extracted from peripheral blood mononuclear cells (PBMCs) using PureLink® Genomic DNA Mini Kit (Invitrogen, USA) according to the manufacturer's protocol. 500 ng of double-stranded DNA was determined by Qubit (Invitrogen, USA) and randomly fragmented to 150-200 bp with Covaris cracker (Covaris, USA). Fragments with specific indexes were hybridized with probes. After PCR amplification and quality control, libraries were sequenced by next-generation sequencing. Agilent liquid phase hybridization was applied to efficiently enrich whole exons which would be sequenced on Illumina platform. Agilent SureSelect Human All ExonV5/V6 (Agilent Technologies, USA) with reagents were used for sequencing libraries and capture, which was recommended by the instruction manual and followed by optimized experimental procedures.

Sequencing was performed on an Illumina HiSeq X sequencer with a paired-end read length of 150 bp in the Genomics Core Facility at Novogene (Genome Sequencing Company, Beijing, China). Data generated in this study will be submitted to the National Center for Biotechnology Information (NCBI) BioProject.

Preparation of a List of Candidate Genes Associated with Rheumatoid Arthritis

A list of 159 candidate RA-associated genetic variants reported by previous genome wide association studies (GWAS) with the P value threshold of P<1×10 −5 , as shown in Table 3, was prepared based on Rheumatoid Arthritis associated genes in the NHGRI GWAS Catalog (Welter D et al., Nucleic acids research 2014; 42:D1001-D1006) and literatures (Freudenberg J et al., Arthritis Rheumatol 2014; 66:1121-1132; Manolio T A et al., Nature 2009; 461:747-753; Okada Y et al., Nature 2014; 506:376-381; and Diogo D et al., The American Journal of Human Genetics 2013; 92:15-27).

Analysis of the Whole-Exome Sequencing (WES)

To analyze the entire cohort of samples for genotype calls, variant analysis and joint genotyping were performed according to the pipeline recommended by the Genome Analysis Toolkit software and the GATK Best Practices procedures on RA patients and healthy controls (San Lucas F A et al., Bioinformatics 2012; 28:421-422; and Dong C et al., Human molecular genetics 2015; 24:2125-2137). Briefly, Burrows-Wheeler Aligner (BWA) software was utilized to align the raw sequencing reads in FASTQ formats to the 1000 Genomes (GRCh37+decoy) human genome reference. The BWA alignment files were converted to BAM files with SAMtools v1.1, which was used for sorting the BAM files. Duplicate reads were marked for BAM files with Picard MarkDuplicates (https://sourceforge.net/projects/picard/). The coverage and depth were computed based on the final BAM file. Local realignment, base quality recalibration, variant calling, joint genotyping, and variant quality score recalibration and filtration were applied using with GATK v3.7 (https://software.broadinstitute.org/gatk/). Default settings were used for BWA, SAMtools, Picard and GATK tools.

Further filtration for the joint genotyped variants was performed using Variant Tools (San Lucas F A et al., Bioinformatics 2012; 28:421-422). The inventors applied the following filters to generate a list of preliminary variants by removing false-positive variants through Variant Quality Score Recalibration with tranche truth sensitivity threshold <99.00, as well as variants with low read depth (DP)<10 and poor genotyping quality (GQ)<20, keeping exonic or splicing variants based on University of California, Santa Cruz (UCSC) genome browser build 37 human Reference Sequence Gene annotation, and removing synonymous variants.

From the preliminary variant list, variants annotated as “pathogenic” in ClinVar and deleterious variants were identified, respectively, including those candidate genes that overlapped with previous studies or passed the case-control gene burden test threshold. Deleterious variants were predicted to be damaging (disease-related, D) or benign/neutral (tolerated, T) based on LR score determined by logistic regression (LR) model (Dong C et al., Human molecular genetics 2015; 24:2125-2137). The novel deleterious variants were divided into the rare and common variant groups, which were distinguished by minor allele frequency (MAF) in Chinese Southern population from the 1000 Genomes Project phase III study.

Analysis of Burden Association Signal

Case-control gene burden analysis was assessed on both rare and common deleterious variants to investigate causal genes using RA patients with >80% Chinese ancestry as cases and two types of controls: 105 southern Chinese samples from the 1000 Genomes Project phase III study and 66 healthy controls with >80% Chinese ancestry. Regardless of DP or GQ, all available genotype calls contributed to the number of allele count across the retained deleterious variants in each individual gene. The gene burden ratio was calculated by dividing the allele frequency in cases by the allele frequency in controls. We identified an enrichment of deleterious variants in a gene according to the gene burden ratio >1.5-fold with both types of controls, or the deleterious alleles in the gene with at least 3 RA cases if zero allele frequency in the controls. The inventors further identified genes with rare variants that were homozygous in RA cases but not present in controls, which were considered greater contribution to functional impact.

Pathway Analysis

To discover enriched functional-related gene groups, pathway analysis was performed using DAVID Bioinformatics Resource 6.8 program (DAVID 6.8) (https://david.ncifcrf.gov/summary.jsp) with a Modified Fisher Exact P value less than 0.05 as the significance threshold and strong enrichment in the annotation categories.

Homology modeling is one of the best and reliable ways to construct the three dimensional (3D) structure of protein (Yamaguchi H et al., European journal of medicinal chemistry 2011; 46:1325-1330). Firstly, protein sequence was imported into the Molecular Operating Environment (MOE) 2015.09 software (Chemical Computing Group Inc., Montreal, Canada) to search an optimal template. The top ranked structure based on the Z score towards the target sequence was selected as the template. Target protein sequence and its corresponding crystal structure coordinates of template were separately loaded and aligned. A series of protein models were independently constructed by using a Boltzmann-weighted randomized procedure (Levitt M., Journal of molecular biology 1992; 226:507-533). Amber force field was applied in the process of construction and energy minimization (Case D et al., Amber 12 reference manual). Finally, the model with the best packing quality function was selected for further full energy minimization, and the stereochemical qualities of protein model was assessed by means of Ramachandran plots.

To analyze the effect on the point mutation in the 3D structure of the protein, the mutant protein were carried out in Residue Scan module of MOE 2015.09 software based on the 3D structure of homology modeling. In addition, we further analyze the hydrogen bonds, solvent interactions, metal ligation and non-bonded interaction between the target mutant residue and its surrounding key amino acid residues.

Results

1. Deleterious Variants in Novel RA Candidate Genes

WES data were generated from 58 RA patients with a median coverage of 76-fold on targeted exome regions ( FIG. 2A ). An average of 96% of all targeted regions was covered by at least 20-fold. The healthy control group had a median coverage of 68-fold on targeted exome regions, and an average of 94% of those regions was covered by at least 20-fold ( FIG. 2B ).

As shown in the flow chart of FIG. 3 , a total of 3,537,952 variants were identified from 58 RA samples. After applying the quality filters and removing synonymous variants, we found 72,024 exonic and splicing variants, including nonsynonymous substitutions and a small number of stop-gain, stop-loss, frameshift and non-frameshift indels (Table 4). Of these, 135 variants were identified as deleterious based on the “pathogenic” annotation in ClinVar, and an additional 21,755 variants predicted to be deleterious were identified using an ensemble logistic regression score.

It was surprising that the identified genes were not found in previously reported candidate risk variants GWAS data (group 1 and group 3 in FIG. 3 ; Table 5 and 6), such as HLA-associated genes. After reviewing the location of the reported RA-associated variants, the inventors found that only 9 of over 200 variants are located in the exome area, including CTLA4 (rs231775), FCGR2A (rs1801274), IL6R (rs2228145), OLIG3 (rs2230926), PTPN22 (rs2476601), RTKN2 (rs3125734), SH2B3 (rs3184504), TNFAIP3 (rs223092) and TYK2 (rs34536443). Since the inventors applied WES technique to focus on exome regions, the results are different from the previously reported studies. It is advantageous that the presently disclosed method of identification of RA associated genes is capable of providing more comprehensive information of RA associated genes. According to the results, a number of new RA associated genes have now been identified and may be useful in developing advanced methods for diagnosis and treatment of RA.

Interestingly, two novel risk variant loci were identified associated with TGFβ1 (transforming growth factor β1) and FOXP3 (forkhead box P3) genes (group 4 and group 6 in FIG. 3 ; Table 5) whereas other known variations of these two genes were previously found to be involved in the risk to RA (Zhou T B et al., J Recept Signal Transduct Res 2014; 34:469-475; and Paradowska-Gorycka A et al., J Rheumatol 2015; 42:170-180).

In order to identify novel genes and pathways that could enhance understanding of RA pathogenesis, the inventors performed a gene burden analysis to identify genes for which deleterious variants were enriched in the Han Chinese RA samples compared to healthy control and public control samples. Six such genes were identified (group 2 in FIG. 3 ; Table 5 and 7). Of these, a missense variant of SAA1 (Serum Amyloid A1) was found in 3 RA patients but not present in healthy controls. SAA1 is highly expressed in response to inflammation and tissue injury, and strongly associated with activity of the disease and risk of cardiovascular and renal involvement in RA patients, suggesting that this novel deleterious variant may potentially contribute to RA disease risk through its interference with pro-inflammatory effectors. Additional pathogenic variant of OXCT1 (3-Oxoacid CoA-Transferase 1) was predicted to be damaging (disease-related, D) in the RA patients, encoding Succinyl-CoA:3-ketoacid coenzyme A transferase 1 (SCOT1), which is a key enzyme for synthesis and degradation of ketone bodies involved in cardiovascular disease.

Rare variants are more likely to predict a significant impact on protein function and result in clinically relevant consequences than common ones. Thus, the inventors grouped variants that were indicated to be deleterious into rare (minor allele frequency <1%) and common variants, which did not overlap with previously reported candidates (Table 5). Performing a gene burden analysis for variants within each of these groups, we identified 241 genes (group 4 in FIG. 3 ; Table 5) with rare, deleterious variants specifically enriched in our Chinese RA samples compared to healthy control and public control samples. Notably, since the functional impact of rare and deleterious variants is likely to be greater when present as homozygote, 5 rare and deleterious homozygous variants (NCR3LG1, RAP1GAP, CHCHD5, HIPK2 and DIAPH2) were identified in the Chinese RA samples and absent in the controls (group 5 in FIG. 3 ; Table 3). Finally, we also identified 144 genes with common and deleterious variants in RA patients (group 6 in FIG. 3 ; Table 5).

2. Pathway Analysis

Using the method as disclosed herein, the inventors identified a total of 381 genes as candidates for increased risk of RA (Table 5). In order to further identify the associated biologic pathways, the inventors performed the functional enrichment analysis using DAVID 6.8 and identified the pathways of the extracellular matrix (ECM)-receptor interaction, protein digestion and absorption, focal adhesion and glycerophospholipid metabolism as significantly overrepresented (Table 8), which were reported to be relevant in pathogenesis of arthritis (Lv W et al., Mol Biosyst 2015; 11:2986-2997; and Choe J Y et al., Rheumatology (Oxford) 2016; 55:928-938).

In order to identify variants that might predispose RA patients to disease duration, the inventors repeated the variant filtration and gene burden analysis on Chinese RA samples with the disease duration 3-year compared to the disease duration ≤1-year. A total of 277 genes were identified (Table 9) compared to the 381 genes identified in the case-control comparison (Table 5). Of these, 87 genes were unique to disease duration with exonic variants (Table 10). Pathway analysis performed on the 87 genes identified olfactory transduction pathway as significantly overrepresented (Table 8), including OR14C36, OR4A15, OR52N4, OR6C74, OR6C75, OR7G3 and OR9K2.

3. Structural Analysis and Function Change Prediction of Potential Biomarkers

In order to gain structure insights of the potential biomarkers with pathogenic variants into the clinical conditions of RA patients, the inventors derived a three-dimensionally structure model of SAA1 Gly90Asp (rs79681911) and SCOT1 Thr58Met (rs75134564) by combining homology modeling with point mutation in MOE 2015.09 package. The crystal structures of human SAA1 protein (PDB code: 41P8.A) and SCOT1 protein (PDB code: 3K6M.C) were selected to be used as templates due to their optimal identity with the target sequences of SAA1 (Protein RefSeq: NP_000322.2) and SCOT1 (Protein RefSeq: NP_000427.1), 83.6% and 83.5%, respectively ( FIG. 4 ). The SAA1 and SCOT1 models with the best packing quality function and full energy minimization were assessed by Ramachandran plots, indicating that the phi and psi backbone dihedral angles in the models were reasonable ( FIG. 5 ).

Structural analysis of SAA1 Gly90Asp ( FIGS. 6A and 6B ) revealed that the substitution of glycine with aspartic acid induced the formation of two pairs of hydrogen bonds with two threonine residues (Ala91 and Asp93), exhibiting more stable structure of loop region and promoting the polar interaction. Moreover, this mutation shortened the length of a helix 4, which may affect the stability of SAA1. Structural changes in SAA1 protein caused the surface of Asp90 to be exposed in solvent environment, leading to the increased hydrophilic region. In addition, the construction of 3D models in SCOT1 and its mutant Thr58Met revealed that this substitution resulted in the disappearance of the hydrogen bond between Thr58 and Asp206, the reduction of intramolecular polar interactions and the expansion of hydrophobic region ( FIGS. 6C and 6D ), suggesting its potential function alteration in RA pathogenesis.

4. Discussion

The inventors performed perspective WES aiming to identify potentially causal biomarkers in a cohort of Chinese RA patients. The inventors used the method as disclosed herein to focus on investigating the occurrence frequency of variants in genes previously associated with RA as well as novel genes. Despite known variants of TGFβ1 and FOXP3 genes associated with increased RA risk, two novel risk variant loci in these two genes were for the first time identified to be implicated in the RA risk (group 4 and group 6 in Table 5). A novel splicing variant (rs199982059) of TGFβ1 was found to be significantly enriched in 4 RA patients, but absent in healthy controls. TGFβ1 is a pivotal protein in the pathogenesis of a number of autoimmune disorders and its dysregulation is also increasingly implicated in the risk of developing RA. RNA splicing is a focal point on connection between genetic variations and complex disorders, and this novel splicing variant of TGFβ1 might provide new insights into the genetic determinants of RA disease. In addition, a novel missense variant (chrX:49114808) of FOXP3 was observed in 8 RA patients. FOXP3 is a unique regulatory T cell (T reg )-specific marker and important in the development of RA-derived T reg cells as a transcriptional factor. In spite of the other known variants in TGFβ1 and FOXP3 genes associated with RA, these two newly-identified variants in our Chinese RA patients may offer the novel genetic contributions to the RA risk.

The inventors have also identified six novel and deleterious genes that are classified as pathogenic in ClinVar database (Table 7). Of these, a missense variant (rs79681911) of SAA1, initially characterized by serum amyloid a variant (OMIM 104750) and required for the amyloidosis disease process, was identified in the RA patients. SAA1 has been reported to play a pathogenic role in the pro-inflammatory cascades in RA, therefore, this novel deleterious variant may be implicated in RA risk as a sensitive indicator of inflammatory activity. Additional pathogenic variant (rs75134564) of OXCT1 was predicted to be disease-related in 4 RA patients based on LR score, which previously implicated in succinyl-CoA acetoacetate transferase deficiency (OMIM 601424) in clinic. OXCT1 encoding enzyme SCOT1 is essential for ketone body metabolism and involved in cardiovascular disease, which are shown to be strongly associated with the course of RA, suggesting this enzyme may potentially contribute to RA prognosis. Importantly, the 3D structural analysis of these two potential biomarkers revealed that the substitution of mutation points may be involved in the functional alteration of the proteins and further impact on RA disease progression ( FIG. 6 ).

The inventors sought to identify novel genes or biological candidate pathways fundamental to the risk of RA disease, including both rare and common variants. To elucidate additive effects of polygenic variants that affect the same gene or pathway, the inventors performed gene burden test and pathway analysis. Notably, the biological impact of rare and deleterious variants is likely to be greater when present as two copies. In the study, 5 homozygous variants (group 5 in FIG. 3 ; Table 5) were detected in the RA patients but not in healthy controls. Intriguingly, a non-frameshift indel variant (rs61406813) of NCR3LG1 (natural cytotoxicity triggering receptor 3 ligand 1) was identified in the RA patients as homozygote. NCR3LG1 could be detected on monocytes and neutrophils after application of inflammatory stimuli, and it was initially described as a tumor cell-expressed ligand of NKp30, which is found to be implicated in RA-associated inflammation. Additionally, a missense variant (rs61014678) of RAP1GAP (RAP1 GTPase Activating Protein) was identified as damaging (disease-related, D) by determinant of LR model in the RA patients. RAP1GAP regulates the activity of the ras-related RAP1 protein, which involves in induction of apoptotic pathway in synovial fibroblasts and plays a critical role in oxidative stress and T cell behavior in RA synovial tissues. Thus, these two homozygous variants may perform stronger functions in RA pathogenic mechanisms.

The WES analysis totally identified 381 genes that may partially contribute to RA pathogenesis and disease progression, including 3 genes (TGFβ1, FOXP3 and SAA1) previously implicated in RA and 378 novel candidate genes. Biologic pathway analysis might help us to deeply understand RA pathogenesis, and previously biological pathways have been identified from genes in large-scale association analysis of GWAS data (Table 6), such as autoimmune thyroid disease, natural killer cell mediated cytotoxicity and T cell receptor signaling pathways. The inventors deciphered enrichment of the identified deleterious genes within additional pathways of ECM-receptor interaction, protein digestion and absorption, focal adhesion and glycerophospholipid metabolism based on our WES data (Table 8), which have been implicated in the autoimmune conditions or pathogenesis of RA. The inventors also sought to identify potential deleterious variants associated with disease duration among RA patients. The pathway analysis focusing on variants enriched among RA patients with disease duration ≥3-year highlighted seven novel genes in olfactory transduction pathway (Table 8), which has been previously reported to be implicated in regulating inflammatory responses.

Pathogenesis of RA is complicated and includes both environmental and genetic factors. Recently, gut microbiota has been evident of being implicated in RA pathogenesis and treatment responses as a critical environmental factor that influences metabolic and immune homeostasis, involvement of protein digestion and absorption, glycerophospholipid metabolism and olfactory transduction pathway, which were also enriched by novel candidate genes identified in the Chinese RA patients (Table 5). In addition, the homozygous variant NCR3LG1 (group 5 in Table 5) may mediate autoimmune and microbial infection-induced inflammation by associating with the ligand of NKp30. Therefore, these involved novel deleterious genes might be convincingly considered genetic contributions to microbial alteration in relation to the pathogenesis and development of RA.

Genetic factors on the X chromosome always contribute to the increased risk of developing autoimmune disorders in females compared with males, such as RA. According to the method as described herein, four novel and deleterious variants were investigated to be associated with sex bias in the Chinese RA patients, including OTC (Ornithine Transcarbamylase) (rs72554348), DIAPH2 (Diaphanous Related Formin 2) (rs363755), ARSE (Arylsulfatase E) (rs56393981) and FOXP3 (chrX:49114808) ( FIG. 7 ). Notably, OTC, ARSE and FOXP3 were previously reported to be implicated in x-linked diseases (Bennett C L et al., Nature genetics 2001; 27:20; Luksan 0 et al., Hum Mutat 2010; 31:E1294-1303; and Jeon G W et al., Ann Clin Lab 2013; 43:70-75), in which these three novel variants identified in the present study are also associated with female, supporting that the association of variants on X chromosome and RA may further provide molecular evidence as a risk factor contributing to increased susceptibility in Chinese female RA patients.

In summary, the inventors have performed WES to present support and improve our understanding of associations with genetic biomarkers that may be involved in the development of RA in the Chinese population. The biomarkers highlighted include previously implicated genes as well as novel genes and pathways, involved in regulation of adaptive immune response, transmission of nerve impulse and chromosome organization. This study significantly extends the work of GWAS and provides new insight into fundamental etiologic mechanisms in this common autoimmune disease. Taken together, these novel biomarkers can be served as novel biomarkers for valid diagnosis tools for identification of RA patients from normal people specifically for Chinese Han population.