Heritability and Genome-Wide Association Analyses of Serum Uric Acid in Middle and Old-Aged Chinese Twins

Serum uric acid (SUA), as the end product of purine metabolism, has proven emerging roles in human disorders. Here based on a sample of 379 middle and old-aged Chinese twin pairs, we aimed to explore the magnitude of genetic impact on SUA variation by performing sex-limitation twin modeling analyses and further detect specific genetic variants related to SUA by conducting a genome-wide association study. Monozygotic (MZ) twin correlation for SUA level (rMZ = 0.56) was larger than for dizygotic (DZ) twin correlation (rDZ = 0.39). The common effects sex-limitation model provided the best fit with additive genetic parameter (A) accounting for 46.3%, common or shared environmental parameter (C) accounting for 26.3% and unique/nonshared environmental parameter (E) accounting for 27.5% for females and 29.9, 33.1, and 37.0% for males, respectively. Although no SUA-related genetic variants reached genome-wide significance level, 25 SNPs were suggestive of association (P < 1 × 10−5). Most of the SNPs were located in an intronic region and detected to have regulatory effects on gene transcription. The cell-type specific enhancer of skeletal muscle was detected which has been reported to implicate SUA. Two promising genetic regions on chromosome 17 around rs2253277 and chromosome 14 around rs11621523 were found. Gene-based analysis found 167 genes nominally associated with SUA level (P < 0.05), including PTGR2, ENTPD5, well-known SLC2A9, etc. Enrichment analysis identified one pathway of transmembrane transport of small molecules and 20 GO gene sets involving in ion transport, transmembrane transporter activity, hydrolase activity acting on acid anhydrides, etc. In conclusion, SUA shows moderate heritability in women and low heritability in men in the Chinese population and genetic variations are significantly involved in functional genes and regulatory domains that mediate SUA level. Our findings provide clues to further elucidate molecular physiology of SUA homeostasis and identify new diagnostic biomarkers and therapeutic targets for hyperuricemia and gout.

The SUA level is mediated by the interplay between genetic and environmental factors. So far, the magnitude of genetic sources of variance in SUA level has been previously explored in several population studies (10)(11)(12)(13)(14)(15)(16). A strong genetic component was indicated with heritability estimates approximately ranging from 35 to 77%. Additionally, the genetic epidemiology has presented an enormous impact on the molecular physiology related to SUA homeostasis by genome-wide association study (GWAS), identifying several genetic loci located in key urate transporters such as SLC22A7, SLC2A9, SLC22A11, SLC22A12, ABCG2, etc., and a number of additional intriguing genetic networks (17)(18)(19).
Although intensively deployed, no GWAS has, to our knowledge, yet been performed on a sample of middle and old-aged Chinese twins. Chinese population differs in the genetic constitutions and a multitude of life style like dietary habit, work type, and physical activity from other ethnic populations in the world. Genetically related individuals, such as twin pairs, would highly confer increased power in genetic association analysis and efficiently identify genetic variants underlying human complex diseases (20).
Based on a sample of 379 middle and old-aged Chinese twin pairs, we explore the magnitude of genetic impact on SUA variation by performing twin modeling analyses and replicate previous findings on heritability of SUA level and further conduct a GWAS to detect specific genetic variants associated with SUA.

Participants
The sample collection was carried out through the Qingdao Twin Registry, and details of study recruitment have been described previously (21,22). Participants who were with gout, systemic lupus erythematosus, eGFR < 60%, or serum creatinine level >1.4 mg/dL were excluded, and incomplete co-twin pairs were also dropped. The final sample consisted of 379 complete twin pairs with a median age of 50 years (95% range: 41-69 years), including 240 monozygotic (MZ) pairs (114 male and 126 female pairs) and 139 dizygotic (DZ) pairs (41 male, 39 female, and 59 opposite-sex pairs).
All co-twin pairs undertook a health examination after a 10-12 h overnight fast and completed a questionnaire. Serum and plasma were separated from blood cells in the field within 30 min and kept frozen at −80°C. The zygosity was determined by using 16 multiple short tandem sequence repeat DNA markers (23,24). SUA level was measured on the Semi-automatic Analyzer (Hitachi 7600, Japan) and transformed following Blom's formula for normality.
This study was approved by the Regional Ethics Committee of the Qingdao CDC Institutional Review Boards. Prior written informed consent was achieved for all participants. The ethical principles of the Helsinki Declaration were followed.

statistical analysis
Heritability Data preparation and descriptive analyses as well as genetic analyses were performed with SPSS version 22.0 and Mx program, 1 respectively. Twin pair phenotypic correlations per zygosity were firstly measured by calculating the Pearson's product-moment correlation coefficients. The higher correlations of MZ than those of DZ twin pairs indicated the genetic effect on individual differences in SUA level.
Then, standard structural equation modeling methods were used for sex-limitation twin modeling based on the classical twin methods. The variation was decomposed into sources of additive genetic (A), common or shared environmental (C), and unique/ nonshared environmental (E) parameters. After the general sex-limitation ACE model was firstly fitted, we then fitted its sub-models: the common effects sex-limitation models by setting the male-specific additive genetic effects (Am ʹ ) and sex-specific common or shared environmental effects (Cm and Cf) to 0 and the scalar sex-limitation models by constraining the variance components for females to be equal to a scalar multiple of the variance components for males, respectively.
In order to choose the best fitting model, the likelihood ratio test was applied to compare the performances between the general sex-limitation model and its sub-models. In the likelihood ratio test, twice the difference in the log likelihoods between models was calculated, and change in chi-square against the change in degrees of freedom were tested. The Akaike's information criterion (AIC) was calculated and a lower AIC indicated a better fit when no statistical difference was observed between two models (26). The covariates of age and body mass index (BMI) were adjusted for in the analysis.

SNPs-Based Genome-Wide Association Study
The association of SUA level with SNP genotypes was tested using the GEMMA (25). The covariates of age, sex, BMI, and the first five principle components were adjusted for in the model fitting. The conventional genome-wide significance level of P < 5 × 10 −8 and suggestive evidence level of P < 1 × 10 −5 for this association were adopted (27). We further conducted functional elaboration of GWAS results and predicted putative causal variants in haplotype blocks, likely cell types of action and candidate target genes of noncoding genome by using online HaploReg v4.1 software 2 (28,29). A set of 25 query SNPs (P < 1 × 10 −5 ) was submitted. The enrichments of cell-type enhancers with uncorrected P < 0.05 were reported.

Gene-Based Analysis
We performed gene-based tests on GWAS summary results by using Versatile Gene-based Association Study-2 (VEGAS2) which uses 1,000 genomes data to model SNP correlations across the autosomes and chromosome X (30,31). In the test, the evidence for association from all SNPs was aggregated within a per gene while correcting for linkage disequilibrium and gene size, and genes showing more signal or strength of association than expected by chance were identified. The SNPs from "1000G East ASIAN Population" were adopted. The P < 2.63 × 10 −6 (0.05/19,001) was considered to be genome-wide significant for the association as 19,001 genes being evaluated.

Gene Sets-Based Analysis
A list of significant genes (P < 0.05) were included to compute the over-represented gene sets in the gene sets-based analysis using the online version of gene set enrichment analysis (GSEA) program 3 (32, 33   1 | Quantile-quantile plot for quality control check and visualizing crude association for genome-wide association study of serum uric acid (SUA) level. The x-axis shows the −log10 of expected P-values of association from chi-square distribution and the y-axis shows the −log10 of P-values from the observed chi-square distribution. The black dots represent the observed data with top hit SNP being colored, and the red line is the expectation under the null hypothesis of no association. Gene at the best SNP is indicated. ( Table 2). The strongest association was detected with rs346750 (P = 2.50 × 10 −7 ) in an intronic region of EXOC3L2 on chromosome 19q13.32. Among these top signals, two chromosomal loci (17q25.3 and 14q24.3) showed nominal association with SUA level as the locus zoom plots illustrated (Figures 3 and 4). On chromosome 17q25.3, seven SNPs (P = 3.25 × 10 −6 -7.29 × 10 −6 ) and two SNPs (P = 1.30 × 10 −6 -8.93 × 10 −6 ) were located at or near TNRC6C and TMC6/TNRC6C-AS1 genes, respectively. At chromosome 14q24.3, four SNPs rs2270073, rs2270074, rs11621523, and rs2159179 (P = 2.77 × 10 −6 -7.59 × 10 −6 ) were positioned within or closest to PTGR2 gene that was involved in terminal inactivation of prostaglandins. The rs2336742 (P = 7.02 × 10 −6 ) was located at the intronic region of ENTPD5 gene, which was involved in the pathway of purine metabolism. The number of SNPs mapping to ZNF410 and FAM161B was five (P = 3.39 × 10 −6 -7.64 × 10 −6 ) and one (P = 7.64 × 10 −6 ), respectively. All the abovementioned genes showed nominal association with SUA level (P < 0.05) from the following VEGAS2 analysis.
As predicted by HaploReg v4.1, two cell-type specific enhancers (uncorrected P < 0.05) of brain angular gyrus (P = 0.003) and skeletal muscle (female) (P = 0.008) were identified for the set of 25 query SNPs (Table S3 in Supplementary Material). Most of the SNPs were located in intronic regions. Several SNPs were detected within regions with promoter histone marks or enhancer histone marks and could change DNA motifs for DNA-binding proteins, and thus would have regulatory effects on gene transcription (Table S4 in Supplementary Material). We compared previously reported 2,368 significant SUA-associated SNPs in a series of studies with our results. Although no genome-wide significant SNPs were identified in our study, we defined our SNPs with P < 0.05 as supportive to the reported SNPs. And 57 SNPs located  in genes SLC2A9, ABCG2, LRRC16A, LOC107986260, GLUT9,  SCGN, LOC107986971, TFCP2L1, TET2, KCNQ1, FRAS1,  SLC16A9, RAF1P1, LOC100129344, LOC107986581, LRRC16A,  LOC100287951, WDR1, PDZK1, and LOC107986260 could be replicated (Table S5 in Supplementary Material).
While no genes achieved genome-wide significance level, a total of 167 genes were observed to be nominally associated with SUA level (P < 0.05) from VEGAS2 analysis. Genes of TNRC6C-AS1, ZNF410, TNRC6C, FAM161B, PTGR2, SESTD1, RGS20, ENTPD5, EXOC3L2, DNAH3, and TMC6 had already been indicated in the SNPs-based analysis ( Table 2), whereas the others were novel. The well-known urate transporter SLC2A9 gene was also identified. The top 20 genes ranked by P-values were listed in Table 3.
In the gene sets-based analysis using GSEA program, one REAC-TOME gene set and 20 GO gene sets (FDR q-value < 0.01) were presented in Table 4. The only REACTOME pathway was transmembrane transport of small molecules (FDR q-value = 0.020). And the GO gene sets were involved in ion transport, regulation of catabolic process, transmembrane transporter activity, hydrolase activity acting on acid anhydrides, cellular process, regulation of immune system process, etc.

DiscUssiOn
In this investigation based on 379 twin pairs, we explored the proportion of genetic sources in SUA variation, and confirmed the genetic variants underlying this trait by GWAS. The MZ twin correlation for SUA level was larger than for DZ twin, indicating the presence of genetic influence (Table S2 in Supplementary Material). In fitting the sex-limitation ACE models, the common effects sex-limitation model of dropping Am ʹ effects (Model II) was favored over the general model (Model I) (AIC = 443.56, P > 0.05), indicating that there was no evidence for sex-specific additive genetic effects. We then considered whether the common or shared environmental effects for males (Cm) (Model III) and further for females (Cf) (Model IV) could also be fixed to 0. However, the goodness-of-fit statistics indicated that these two models provided a significantly worse fit (P < 0.05). Thus, these two models were rejected and model II remained the favored one. Finally, we considered the scalar sex-limitation models (Model V and Model VI) by constraining the variance components of females to be equal to scalar multiples of the variance components of males. Both of the models provided a significantly worse fit than Model II (P < 0.05). Hence, we concluded that the Model II was the best fitting model, in which additive genetic parameter (A) explained larger proportion of SUA variation for females (46.3%) than for males (29.9%), whereas the environmental parameters together (C and E) explained smaller proportion (53.7 vs. 70.1%) ( Table 1).
The sex-difference in the genetic and environmental effects on SUA variation obtained here was in line with the previous Boyle et al. 's genetic study. Based on a sample of 112 twin pairs, they also found a more significant genetic component in control of SUA variation in females than males, whereas a stronger role of environmental component in males (34). We speculated that the genetic architecture may indeed differ across sexes because of the sex differences in selective pressures during human evolution. Additional data from adopted twins and siblings reared together may be used to explore this hypothesis further.
Although no genome-wide significant SNPs were identified in GWAS, we found two promising genetic regions on chromosome 17 around rs2253277 (Figure 3) and chromosome 14 around rs11621523 (Figure 4). The PTGR2 and ENTPD5 genes around the rs11621523 have been emphasized for their roles in SUA level. Prostaglandin reductase 2 (PTGR2) is the enzyme involved in terminal inactivation of prostaglandins (35) which   Catalysis of the transfer of a solute from one side of a membrane to the other, up its concentration gradient. The transporter binds the solute and undergoes a series of conformational changes. Transport works equally well in either direction and is driven by a chemiosmotic source of energy, not direct ATP coupling. Chemiosmotic sources of energy include uniport, symport, or antiport may contribute to renal uric acid metabolism (36,37). Four highly correlated SNPs showing suggestive evidence of association with SUA level were detected within or near PTGR2 gene. The rs2270073 and rs2270074 were detected within a region with promoter histone marks in the 5'-UTR of PTGR2 gene and could change DNA motifs for DNA-binding proteins, which provided strong evidence for their regulatory effects on gene transcription (Table S4 in Supplementary Material). Additionally, rs11621523 and rs2159179, which were located in an intergenic region at 14q24.3 and closest to PTGR2 gene, should also be candidates to be further studied. The protein encoding by ENTPD5 gene was involved in the pathway of purine metabolism in which SUA was a by-product of oxidation. As SNP rs2336742 was located at the intronic region of ENTPD5 gene and could change its DNA motifs, it might be associated with purine metabolism as well as SUA level (Table S4 in Supplementary Material). However, the association of novel TNRC6C, TMC6, and TNRC6C-AS1 genes around the other SNP rs225327 with SUA level still needs to be validated. Finally, the enhancer of skeletal muscle was predicted by submitting the set of 25 query SNPs to HaploReg v4.1 (Table  S3 in Supplementary Material). And the relationship between SUA level and skeletal muscle strength/volume has been fully researched currently (38)(39)(40)(41).
As additional replication, we compared the SUA-associated SNPs reported in a series of studies with ours (18, 42-59) ( Table   Table 4 | Continued S5 in Supplementary Material). A list of SNPs could be replicated, especially the variants located in the well-known SUA-associated genes SLC2A9 and ABCG2. Notably, two SNPs rs2231142 (ABCG2) and rs10008015 (TET2) were also found by Yang et al. (58) and three SNPs rs11996526 (LOC107986971), rs4848700 (TFCP2L1), and rs179785 (KCNQ1) by Li et al. (49) in Chinese populations.
Even though none genome-wide significant genes were found, a total of 167 genes were observed to be nominally associated with SUA level (P < 0.05) from VEGAS2 analysis. The results of GO gene sets from using GSEA program indicated that these genes might be associated with process of generation, catabolism, transport, intracellular signal transduction, and hydrolase activity during SUA metabolism. Besides, several genes were enriched in pathway of transmembrane transport of small molecules, including solute carrier family (SLC14A2, SLC22A1, SLC2A9, SLC5A3, SLCO1B3, and SLC12A5) and ATPase family (ATP6V1H and ATP11B), which strengthened their significance in regulating SUA transport process and thus further influencing SUA level ( Table 4). Except for the PTGR2 and ENTPD5 genes being abovementioned, the GPR151 gene should be also noted. This gene encodes an orphan member of the class A rhodopsinlike family of G protein-coupled receptors (GPCRs) and thus influences the GPCRs activity. And the GPCRs could regulate the assembly of a multienzyme complex for purine biosynthesis (60).
The association of well-known urate transporter gene SLC2A9 (P = 0.001) with SUA level has previously been reported (61,62). Other genes were of unknown function in terms of SUA level or purine metabolism currently, whereas they may also be interesting potential candidates to be future researched and validated, especially the top 20 genes ( Table 3).
Several strengths must be noticed in our study. First, our results based on the twin data of SUA level would be credible. Phenotype variation may be under the effect of subjects' genetic background, age, gender, and environmental exposures as well as some experimental variables related to sampling, processing, and data analysis. Genetically related individuals, such as twin pairs, would highly confer increased power in genetic association analysis and efficiently identify genetic variants underlying human complex diseases (20). Second, given the various genetic constitutions and multitude of life style among different ethnic populations in the world, this is the first GWAS conducted in the sample of middle and old-aged Chinese twins.
Nevertheless, our study has potential limitations as well. First, our study was with relatively small sample size and limited statistical power resulting from the challenges of identifying and recruiting qualified twin pairs. The results presented here, however, provided a useful reference for hypotheses to be further replicated and validated for exploring increased SUA level. Given the genetic effect on SUA variation is expected to comprise a large number of SNPs possessing very small effect size, a meta-analysis with larger samples will be desirable and ideal. Second, even though we replicated parts of our GWAS results by comparing with results generated from external and independent datasets, most of the SNPs didnot reach the genome-wide significance level. In addition, as no other study has explored the differential expressed genes based on SUA-discordant samples, we cannot yet validate our findings further.
In summary, we have confirmed that genetic factors are significant in explaining SUA level variability through twin modeling. Two novel suggestive regions located at chromosomes 17 and 14 were identified. Twenty-five SNPs reached suggestive evidence level of association with SUA and most of them could have regulatory effects on gene transcription, and 167 genes nominally associated with SUA level were involved in significant biological functions related to uric acid generating and metabolism. The potential candidate biomarkers of SUA level reported here should merit further verifications.

DaTa aVailabiliTY sTaTeMenTs
The SNPs datasets for this study have been deposited in the European Variation Archive (EVA) (Accession No. PRJEB23749).

eThics sTaTeMenT
All participants provided written informed consent for participating in the study which was approved by the local ethics committee at Qingdao CDC, Qingdao, China.
aUThOr cOnTribUTiOns DZ and WW contributed to the conception and design. HD and CX organized the collection of samples and phenotypes. WW and YW contributed to sample data and sequencing data management. QT and SL analyzed the sequencing data and WW and CX interpreted the analysis results. WW and DZ drafted the manuscript, YW and HD were involved in the discussion, and QT, SL, and DZ revised it. All the authors read the manuscript and gave the final approval of the version to be published. All the authors agreed to be accountable for all aspects of the work.

acKnOWleDgMenTs
We thank Weilong Li and Jesper Lund for their help on the installation and application of related software and the corresponding results interpretation.

FUnDing
This study was supported by the grants from the National Natural Science Foundation of China (31371024), the EFSD/CDS/Lilly Programme award (2013), and the project of training and studying abroad for full-time doctoral candidate in Faculty of Medicine in Qingdao University.