Patients (n=10) who were attending a clinic for arthritis were recruited between July 2016 and June 2017. Additionally, three families were recruited in the present study as: i) Presenting with autosomal dominant inheritance and ii) healthy family members agreed to participation in the present study. Patient information is summarised in Table I. The present study was conducted according to the principles outlined in the Declaration of Helsinki. The study protocol was approved by the Ethics Committee of The First Affiliated Hospital, Wenzhou Medical University (Wenzhou, China; approval no. 2018-020). Written informed consent was obtained from all participants or their guardians. Patients with gout were clinically evaluated by physicians based on the 2015 gout classification criteria by the American College of Rheumatology/European League Against Rheumatism Collaborative Initiative (21). Patients with inherited metabolic disorders, including Lesch-Nyhan syndrome, were excluded from the present study. Autoantibodies and human leukocyte antigen B27 (HLA-B27) tests were negative in all patients and X-ray analyses of the affected joints were performed to exclude other potential diseases. The classical symptoms of patients with gout are characterized by the rapid development of monoarticular arthritis, which is accompanied by swelling and redness of the first metatarsophalangeal joint (MTP1).

The human genome hg19 was used as the reference genome. Genomic DNA was extracted from peripheral blood leukocytes. In total, 2 µg genomic DNA from each sample was sheared to fragment with a length of 150 base pairs (bp) using the Covaris S220. Subsequently, ligation of small fragments of DNA with A-tail and adapters was conducted. A genomic library was constructed subsequent to the amplification of adapter-ligated DNA using an Agilent SureSelect Library Prep kit (Agilent Technologies, Inc.) according to the manufacturer's protocol, and the samples from each individual were marked with a unique index. Whole-exome capture was performed using the Agilent SureSelect Human All Exon v5 kit (Agilent Technologies, Inc.) according to the manufacturer's protocol. High throughput sequencing was performed using an Illumina HiSeq 4000 sequencer (Illumina, Inc.).

All raw sequencing data obtained from these three families were processed in a similar manner, according to a customized bioinformatics pipeline (22). To remove sequence adapters and low-quality reads, the raw reads were filtered using the FastQC software program, version 1.11.4 (http://www.bioin-formatics.babraham.ac.uk/projects/fastqc/). The filtration criteria excluded Phred-scaled quality scores <30 and read lengths <80 bp. Subsequent to removing the low-quality reads, the remaining reads were used for further analysis. FastQ reads were aligned to the human reference genome (GRCh37/hg19) using the Burrows-Wheeler alignment tool (23) and further visualized using SplicingViewer software (24). Then, Genome Analysis ToolKit (GATK; version 4.0.10.0; https://gatk.broadinstitute.org/hc/en-us) was used to remove duplicated reads and reads mapped to multiple genome locations. In addition, local realignment and map quality score recalibration were performed. Candidate variants were then identified using the GATK Unified Genotype (version 4.0.10.0; https://gatk.broadinstitute.org/hc/en-us).

mirTrios with an integrated ANNOVAR tool were used to annotate all the detected mutations according to an in-house pipeline (25). The minor allele frequency (MAF) was obtained for each variant from publicly available databases (25), including ExAC, UK10K, dbSNP147, 1000 Genomes and ESP6500, and from in-house exome data. The detected variants with a MAF >0.01% present in any of the aforementioned databases were removed (26). Subsequently, the effects of the detected variants were predicted using four tools: i) Polymorphism Phenotyping v2 (PolyPhen2) (27); ii) Likelihood Ratio Test (LRT) (https://evomics.org/resources/likelihood-ratio-test/); iii) MutationTaster (28); and iv) Functional Analysis through Hidden Markov Models (FATHMM) (29). A missense variant was considered deleterious if the variant was predicted to be deleterious or damaging by ≥3 of the four genetic prediction tools. The remaining variants were considered to be high-confidence causative variants.

The amino acid sequence of human protein kinase CGMP-dependent 2 (PRKG2) was retrieved from Uniprot (http://www.uniprot.org/uniprot/Q13237). Crystal structures of the wild-type PRKG2 protein were obtained from the protein data bank (PDB; http://www.rcsb.org/structure/5BV6) and visualized using PyMol software (version 1.8) (https://pymol.org/2/). Homology modelling of mutated PRKG2 protein structures was performed using SWISS-MODEL (https://swissmodelexpasy.org/.) (30). In total, 120 residues were modelled in the PRKG2 structure, including 102 residues in the cGMP-binding region.

Various candidate genes for hyperuricemia and gout were reported in previous studies. The literature search was performed in June 2017 by searching 'gout' AND 'genes' in PubMed (https://www.ncbi.nlm.nih.gov/pubmed). Then all genes in all literature obtained were included. Following the literature search, the candidate genes associated with hyperuricemia and gout were collected and are listed in Table SI.

Temporally rich transcriptome data extracted from the Genotype-Tissue Expression (GTEx) project (https://gtexportal.org/home/, accessed January 2018) were used to build the co-expression network. The Pearson correlation coefficients (r) for the gene co-expression levels were calculated for each pairwise combination between different genes. To investigate similarities among the genes forming the PPI network and analyse their functions, significantly enriched DO terms were identified using the R package DOSE (version 2.0) and a hypergeometric test (31). To assess function similarities between previously reported genes (Table SI) associated with gout and the core risk genes identified in the present study, a PPI network was built. PPI data downloaded from STRING V10 (https://string-db.org) (32) were used for PPI network analysis. To further investigate the gene functions in the PPI network, biological processes (BP) analysis was conducted using ClueGO v2.3.3 (33), a plug-in of Cytoscape.

The patient cohort for the present study was composed of three families (Figs. 1A, 2A and 3). All patients experienced acute monoarticular arthritis affecting the MTP1 and/or knee. The affected joint made walking difficult, was painful to the touch and was occasionally accompanied by fever. Tophus was observed in some of the patients (Fig. 1B), and metabolic disorders, including hypertension, diabetes and hyperlipemia, were identified. The clinical data are presented in Table I. The serum urate level was >480 µmol/l in most cases, and in numerous cases it was >600 µmol/l. Synovial fluid from a number of the patients presented MSU crystals and erosions based on conventional radiography of the affected joint. Autoantibodies and HLA-B27 tests were negative in all patients (Table SII). All affected patients were diagnosed with gout.

To determine the genetic etiology of these families, WES was performed in the three probands. In total, ~13.53 Gb high-quality data was obtained on mean subsequent to removing sequencing adapters and low-quality sequences (Table SIII). For each sample, ≥98.68% of the high-quality data was matched with the human reference genome Hg19. Effective sequence on target was >4.91 Gb, with a minimum of 39.10% fraction of effective bases on target following the removal of polymerase chain reaction duplications. The mean sequencing depth for each sample was 110.16-fold, with >99.00% of target regions being covered at a 4-fold sequencing depth, 98.50% at 10-fold depth and 97.20% at 20-fold depth. Collectively, the quality control data demonstrated a high reliability, which was fundamental for the subsequent analyses.

Subsequent to removing low-quality reads, adapters and duplicated reads from the raw sequencing data, a total of 541,954 single nucleotide variations (SNVs) and 84,415 indels were identified using the GATK tool, which included 67,707 SNVs and 2,200 indels in exonic and splicing regions. Subsequent to applying variant filtration against multiple databases, the number of rare SNVs and indels causing protein change with MAF <0.001 was reduced to 228 and 29, respectively. As a result, following effect prediction by Polyphen2, LRT, Mutation Taster and FATHMM, the variants predicted to be deleterious by >2 prediction tools were validated by Sanger sequencing. Finally, three SNVs in the coding regions of three gout-associated genes were retained and confirmed following Sanger sequencing validation (Table II).

Analysis of WES data indicated a missense mutation (p.Lys157Glu) in the ABCG2 gene in the proband F1:II:1 of family 1 presenting typical symptoms of gout (Fig. 1 and Table I). The ABCG2 protein is involved in the excretion of urate from the intestine and kidney, and its dysfunction causes extrarenal urate underexcretion type (34) and/or renal urate underexcretion type gout (35). The heterozygous missense mutation (c.469A>G) causing a lysine (Lys) to glutamic acid (Glu) substitution was located at the amino acid 157. Subsequently, Sanger sequencing confirmed the presence of this mutation in the patient's affected father while his unaffected mother and sister did not present the A to G change at cDNA nucleotide 469 (Fig. 1A). The MAF of p.Lys157Glu was 8.24×10⁻⁶ and was predicted to be deleterious by all four effect prediction tools (Table II). Additionally, the residue 157Lys is highly conserved among different vertebrate species (Fig. 1C). Subsequently, the tissue specificity of the expression of the human ABCG2 was examined in all major tissues and organs using the Human Protein Atlas database (https://www.proteinatlas.org). ABCG2 exhibited the most abundant protein expression in the small intestine (Fig. 1D). The results reflected abundant mRNA expression in the luminal membrane of the intestine with 136.8 transcripts per million (TPM).

Another heterozygous missense mutation (c.752A>G; p.Asn251Ser) located in exon 5 of PRKG2 gene was detected in proband F2:II:1 and his affected brother F2:II:9, as well as in the patient III9 (Fig. 2A). cGMP-dependent protein kinase 2 (cGKII)/PRKG2 is involved in the regulation of aldosterone and renin secretion (36). In total, nine unaffected family members presented wild-type alleles (Fig. 2A). Importantly, this mutation was not previously observed in any public database. The residue 251Asn was revealed to be evolutionarily conserved (Fig. 2B) and PRKG2 expression was identified to be enriched in the small intestine (18.36 TPM), mainly in glandular cells (Fig. 2C). In the F3, the proband F3:II:1 and his two affected sons shared the same missense mutation consisting of a G to C substitution (c.12G>C) in ADRB3 (Fig. 3). ADRB3 is part of the adrenergic system, which is involved in the regulation of lipid metabolism and glucose homeostasis (37). This single-nucleotide change resulted in a non-synonymous substitution (p.Trp4Cys). PolyPhen2 predicted that this variant was likely damaging (Table II). Sanger sequencing validated that six healthy family members presented homozygous wild-type alleles for ADRB3 (Fig. 3). Similarly to the aforementioned two residues, this position was predicted to be evolutionary conserved by the GERP++ tool (https://omictools.com/gerp-tool).

Previous GWAS studies identified single nucleotide polymorphisms (SNPs) in these three genes that were associated with gout (38-40). All SNPs and three candidate mutations identified in the present study were mapped to schematic representations of ABCG2, PRKG2 and ADRB2 genes (Fig. 4A). Except for the two mutations in PRKG2, all the other mutations were located in the protein functional domains. Furthermore, p.Lys157Glu mutation in ABCG2 was located in the same highly conserved ABC transporter domain as the pathogenic missense variant Gln141Lys, which has been revealed to be associated with an increase in serum uric acid levels (41). The Asn251Ser mutation was in the cGMP-binding region of PRKG2; however, rs768867227 and rs10033237 SNPs, which were previously reported to cause gout, were present in non-coding regions (42).

In order to further study the functional defects caused by mutations in the protein structure, the potential structural differences between the wild-type and mutant proteins were investigated. Therefore, numerous differences were identified between the wild-type and mutant PRKG2 proteins. The wild-type structure (5BV6) was downloaded from the PDB. In the wide-type protein, there were seven hydrogen bonds at residue 251 (Asn); one bond was revealed between Asn251 and Glu191, and the other six connected Asn251 to water molecules (Fig. 4B). The mutated structural modelling revealed the formation of two different hydrogen bonds, one between the mutated Ser251 and Thr250, and another between the mutated Ser251 and Tyr189 (Fig. 4B).

Previous studies have identified that gout-associated genes are closely associated with urate excretion (10-12). To further investigate the expression pattern of the three candidate genes (ABCG2, PRKG2 and ADRB2) identified in the present study and to examine whether their expression patterns were similar to other gout-associated genes identified in previous studies, the transcriptomic data in the small intestinal tissue from the GTEx project were analysed. Based on the co-expression network, these three genes demonstrated similar expression patterns with numerous gout-associated genes, with a Pearson correlation coefficient ranging between 0.602 and 0.898 (Fig. 5A). Among these three genes, the number of candidate genes co-expressed with the ABCG2 gene was the largest (n≤39). Furthermore, to identify functional similarities among these co-expressed genes that may have a function in the development of gout, significantly enriched DO terms were identified using the R package DOSE. Among the top ten most statistically significant DO terms with Bonferroni corrected P<0.05 (Fig. 5B; Table SIV), the majority of the terms were associated with cardiovascular and metabolic diseases, supporting the idea that gout is often accompanied by metabolic syndrome. The two directly associated terms, hyperuricemia [false discovery rate (FDR)=8.02×10⁻¹⁴, hypergeometric test] and gout (FDR=2.32×10⁻¹², hypergeometric test) were the most statistically significant, indicating the significant enrichment of the co-expressed genes in hyperuricemia and gout.

To investigate the association between the three candidate genes and other gout candidate genes, data from the STRING v10 database were analysed. The STRING database contains PPIs, including physical and functional associations. As a result, a total of 79 candidate genes were included in the PPI network (Fig. 6A) with the interaction score ranging between 150 and 967. Solute carrier family 22 member 1 exhibited the strongest interaction with ABCG2, which is expressed in the kidney and mediates the transport of xenobiotics, endogenous organic anions and urate (43). PRKG2 exhibited the strongest interaction with inositol 1,4,5-trisphosphate receptor type 1, which was previously reported to be associated with high serum uric acid concentrations (44). In the next step, all 79 candidate genes in the PPI network were used in the BP with ClueGO plugin. A total of 40 enriched BP terms were divided into eight groups (Fig. 6B) and terms in the same group had similar biological functions. In total, two groups served direct functions in the metabolism of uric acid, and the terms in the two groups with the most significant statistical significance were 'purine nucleoside monophosphate biosynthetic process' and 'urate transport'. 'Glycometabolism' and 'anion transport' groups were also identified following data enrichment analysis. The detailed results of the BP analysis are presented in Table SV.