Whole-exome sequencing and genome-wide association studies identify novel sarcopenia risk genes in Han Chinese.

Sarcopenia is a complex polygenic disease, and its molecular mechanism is still unclear. Whole lean body mass (WLBM) is a heritable trait predicting sarcopenia. To identify genomic loci underlying, we performed a whole-exome sequencing (WES) of WLBM variation with high sequencing depth (more than 40*) in 101 Chinese subjects. We then replicated in the major findings in the large-scale UK Biobank (UKB) cohort (N = 217,822) for WLBM. The results of four single-nucleotide polymorphisms (SNPs) were significant both in the discovery stage and replication stage: SNP rs740681 (discovery p = 1.66 × 10-6 , replication p = .05), rs2272303 (discovery p = 3.20 × 10-4 , replication p = 3.10 × 10-4 ), rs11170413 (discovery p = 3.99 × 10-4 , replication p = 2.90 × 10-4 ), and rs2272302 (discovery p = 9.13 × 10-4 , replication p = 3.10 × 10-4 ). We combined p values of the significant SNPs. Functional annotations highlighted two candidate genes, including FZR1 and SOAT2, that may exert pleiotropic effects to the development of body mass. Our findings provide useful insights that further enhance our understanding of genetic interplay in sarcopenia.

INTRODUCTION

Sarcopenia was defined as age‐related loss of skeletal muscle mass and function, which was associated with a range of other diseases and health problems such as osteoporosis, obesity, and increased mortality (Brzeszczyńska et al., 2017; Cederholm & Morley, 2015; Chen et al., 2014; Fielding et al., 2011; Hayashi, Abe, Fujita, Okai, & Ohira, 2018; Joglekar, Nau, & Mezhir, 2015; Monaco, Castiglioni, Monaco, & Tappero, 2017; Valter, Andrea, Massimiliano, & Marco, 2014); 13%–24% of the elderly white people under 70 years old have sarcopenia, and prevalence is high as 50% or more among those aged 80 years old or older (Baumgartner et al., 1998). The decline in the contractility of lean mass also leads to a decrease in the load on bones, which makes the bones in a state of disuse for a long time, and is easy to induce osteoporosis (Buckwalter, 1995). The diagnosis of sarcopenia is based on an assessment of skeletal muscle mass and/or muscle strength, as well as muscle functions (Zembroń‐Łacny, Dziubek, Rogowski, Skorupka, & Dąbrowska, 2014). Skeletal muscle is under strong genetic control, with heritability estimates of 30%–85% for muscle strength and 50%–80% for muscle mass (Arden & Spector, 2010; Thomis et al., 2004). Whole lean body mass (WLBM) is composed of skeletal muscle (60%), viscera, and some other connective tissues (Vandewoude, Alish, Sauer, & Hegazi, 2012; Wolfe, 2006), with heritability over 50% (Hsu et al., 2012), which is a heritable trait predicting sarcopenia. WLBM can be measured accurately by dual‐energy X‐ray absorptiometry (DXA). Compared to other measurements, DXA is a simple, easy‐to‐use, approachable method with a low dose of radiation and a high degree of precision and accuracy, especially in determining the composition of the body at all local and random sites.

Bone and muscle have persistent paracrine crosstalk throughout life, and they both respond to paracrine and endocrine stimuli in common. A decrease in muscle contraction leads to a decrease in bone load, which leads to a loss of skeletal integrity. Osteoblasts and muscle cells share a common mesenchymal precursor during embryonic development (Hsu et al., 2012). Based on the biological information of the interaction between bone mineral density (BMD) and LM, it is predicted that there may be shared genetic background underlying both traits, that is, pleiotropic genes. Besides, sarcopenia and obesity interact and influence each other on physiological and behavioral levels (Roubenoff, 2000). According to the inconsistent definitions of sarcopenic obesity and studies with different populations, the prevalence of sarcopenic obesity was estimated to be 4%–12% (Cauley, (2015). Both sarcopenia and obesity have high genetic susceptibility, so there may be pleiotropic loci to influence the risks of both diseases.

In the present study, aiming to achieve a better understanding of the genetic etiology of sarcopenia, we conducted a genome‐wide association study of WLBM using the 101 Chinese subjects as discovery sample and a follow‐up replication study in the UK Biobank (UKB) sample. In addition, we conducted a serious of bioinformatic analysis including cis‐eQTL analysis, to explore the functional relevance of the identified genes including FZR1 (OMIM 603,619) and SOAT2 (OMIM 601,311).

MATERIALS AND METHODS

Subjects

The discovery sample consists of 101 unrelated Chinese adults, collected from more than 2000 random subjects in the form of extreme sampling, including 51 low hip BMD individuals (BMD < 0.785 g/cm2）and 50 high hip BMD individuals (Hip BMD > 0.98 g/cm2). All subjects were Han Chinese residents living in the vicinity of the central and western China, Changsha, and Xi'an, and their neighboring areas. All individuals with chronic diseases of vital organs, severe endocrine diseases, severe metabolic diseases, and severe nutritional diseases are not included in the research project. Individuals whose long‐term medication may affect bone mass, bone structure, or metabolism or muscle development were not included in the study.

We replicated the pleiotropic single‐nucleotide polymorphisms (SNPs) identified in the present study by the UKB cohort. In brief, the UKB cohort is a large prospective cohort study of ~500,000 participants from across the United Kingdom, aged between 40 and 69 at recruitment. Ethics approval for the UKB study was obtained from the North West Centre for Research Ethics Committee (11/NW/0382), and informed consent was obtained from all participants. This study used the data requested under the UKB application number 41,542, which was covered by the general ethical approval for the UKB study.

Phenotype measurements and modeling

All subjects completed a structured questionnaire including lifestyle, medical history, family information, anthropometric variables, etc.

In the discovery sample, WLBM and fat body mass (FBM) were measured with a Hologic QDR 4,500 W DXA scanner (Hologic Inc.). Body weight was measured in light clothing on a calibrated balance beam scale. Height was obtained using a calibrated stadiometer. Covariates (including gender, age, age2, height, height2, FBM, and the first five principal components derived from genome‐wide genotype data) were screened for significance with the step‐wise linear regression model implemented in the SPSS software. The WLBM was adjusted through significant covariates, and the residuals were normalized by inverse quantiles of standard normal distribution.

In the replication sample, body composition was quantified by bioelectrical impedance approach. WLBM was approximated by the sum of fat‐free mass at whole body. Phenotype modeling of WLBM was similar with that in the discovery samples. The difference is the phenotype in UKB was mandatorily adjusted by the top 10 principal components to adjust the potential population structure. Association was again examined by the linear regression model.

Genotyping

In whole‐exome sequencing (WES), the qualified genomic DNA subjects were randomly fragmented by ultrasonoscope (Covaris S2). The size of the library fragments is mainly distributed between 250 and 300 bp. Adapters were ligated to both ends of the resulting fragments. After extraction, we amplified DNA by ligation‐mediated PCR (LM‐PCR), purified, and hybridized to the NimbleGen 2.1M human exome array for enrichment. Non‐hybridized fragments were then washed out. Both non‐captured and captured LM‐PCR products were subjected to quantitative PCR to estimate the magnitude of enrichment. Each captured library was then loaded onto Hiseq2000 platform. High‐throughput sequencing was performed for each captured library to ensure that each subjects met the desired average sequencing depth (40× coverage). Raw image files were processed by Illumina base‐calling software 1.7 for base‐calling with default parameters. The sequences of each individual were generated as 90 bp pair‐end reads. With the criteria that quality score of consensus genotype is at least 20, about 20 M reads passed the quality control (QC) in each subject.

Genome‐wide genotypes for all subjects were available at 784,256 genotyped autosome markers, and were imputed into UK10K haplotype, 1,000 Genomes project phase 3, and Haplotype Reference Consortium (HRC) reference panels. All the included subjects are those who self‐reported as white (data field 21,000). Subjects who had a self‐reported gender inconsistent with the genetic gender, who were genotyped but not imputed, or who withdraw their consents were removed. A set of unrelated subjects were then sampled for subsequent analysis with KING (Manichaikul et al., 2010), based on the genome‐wide genotyped variants.

Genotype quality control

We followed strict genotype QC procedure at both individual and SNP levels. At the individual level, sex compatibility was checked by imputing sex from X‐chromosome genotype data with PLINK (Purcell et al., 2007); individuals of unclear sex or of inconsistent with reported sex were removed. At the SNP level, the SNP that violates the Hardy–Weinberg equilibrium (p value < 1.0 × 10–5) was removed. At the same time, the SNP with the Mendel error is set to the missing value.

Population stratification

To monitor and correct for potential population stratification effects, we adopted a principal component analysis (PCA)‐based correction for the detection of the discovery sample. PCA was applied to each individual sample (Zhang, Li, Pei, Liu, & Deng, 2009), and no population outliers were observed.

Meta‐analysis

Fisher's method, also known as Fisher's combined probability test, is a technique for data fusion or "meta‐analysis" (analysis of analyses). It was developed by and named for Ronald Fisher. Significant genes found in the discovery subjects were further replicated in the replication subjects. Fisher's method was then used to analyze two gene‐based association signals jointly. Specifically, Fisher's method was statistically calculated as: where p 1 and p 2 refer to two gene‐level p‐values. Under the null hypothesis of no association, this statistic approximately follows the Chi‐square distribution with 4 degrees of freedom. Note that, Fisher's method is always valid regardless of whether the directions of two effect sizes are consistent.

Functional annotation

Functional annotation of the discovered SNPs was performed using the bioinformatics software HaploReg (Ward & Kellis, 2012). The HaploReg database provides functional information for non‐coded SNPs with multiple functional areas. The functional categories include conservation sites, DNase hypersensitivity sites (DHS), transcription factor binding sites (TFBS), promoter sites, enhancer sites, and others. We annotated these significant SNPs and their neighbor SNPs with strong linkage disequilibrium (LD) pattern (r 2 > .8).

To investigate the association between the identified SNP polymorphisms and the nearby gene expressions, we used the Genotype‐Tissue Expression (GTEx) Portal project dataset (Lonsdale et al., 2013; Westra et al., 2013a) and the Westra et al.'s study (Westra et al., 2013a) to perform cis‐eQTL analysis. The GTEx project collected and RNA‐sequenced multiple human tissues (up to 11,614) from donors who were also densely genotyped, and analyzed associations between SNPs and global RNA expression within individual tissues. We downloaded the summary results of skeletal muscle tissue from the GTEx website (V7) (https://www.gtexportal.org/home/). The Westra et al.'s study performed the largest expression quantitative trait locus (eQTL) meta‐analysis so far reported in non‐transformed peripheral blood samples of 5,311 individuals, with replication in 2,775 individuals (Westra et al., 2013b). We obtained cis‐eQTL results from the research website (http://www.genenetwork.nl/bloodeqtlbrowser/). Cis‐regulation refers to a potential functional genomic variation within or near a gene that significantly affects gene expression of a particular gene.

For candidate genes, we annotated them by constructing gene interaction networks with STRING (Szklarczyk et al., 2015). STRING uses information based on gene co‐expression, text‐mining, and others, to construct gene interactive networks.

Mouse model survey

We surveyed mouse knockout models to evaluate muscle‐related phenotypic consequence of the identified candidate genes, including changes in LBM (p < .05). We primarily searched the International Mouse Phenotyping Consortium (IMPC) database via its web portal (http://www.mousephenotype.org/) (Dickinson et al., 2016).

RESULTS

The basic characteristics of the subjects used in discovery sample are summarized in Table 1.

Table 1

Basic characteristics of the study subjects

	Discovery sample (101 Chinese)
	Male	Female
No. of subjects	50	51
Age	29.32 (3.56)	28.80 (4.70)
Height (cm)	168.31 (5.38)	158.16 (5.89)
Weight (kg)	62.18 (9.78)	52.74 (7.54)
Fat body mass (kg)	16.10(5.20)	26.94(5.37)
Lean body mass (kg)	51.26 (6.96)	38.03 (4.79)

WES generates 332,246 SNPs that are included in the SNP‐based association test in the discovery samples. The Bonferroni correction is used to declare the genome‐wide significance level (GWS, 0.05/332,246 = 1.5 × 10–7). No SNPs achieved GWS. The most significant SNP is rs740681 (p = 1.66 × 10–6). rs740681 is at 19p13.3; previous studies have reported that 19p13.3 was associated with LBM rs7249081 (2 × 10–14) and rs546547223 (5 × 10–6) (Hübel et al., 2019; Tachmazidou et al., 2017).

After adjusting phenotypes by PCA in each individual study, the genomic control inflation factor of the meta‐analysis was 1.142, implying limited effect of potential population stratification. A logarithmic quantile–quantile (QQ) plot of the adjusted test statistics shows a marked deviation in the tail of the distribution, implying the possible existence of true associations (Figure 1).

Figure 1

QQ plot. Logarithmic quantile–quantile (QQ) plot of the discovery sample

We replicated the suggestive SNPs (p < 1.0 × 10–3) in the UKB sample. Four SNPs were significant (p ≤ .05) in the UKB sample (Table 2). The effect direction at the SNP is consistent with that in the discovery stage for phenotypes. The main association results of the consistent SNP are listed in Table 2.

Table 2

Significant association results for SNPs

						Discovery (N up to 101)			Replication (N up to 217,822)			P combined
						Beta	SE	P	Beta	SE	P
SNP	Chr	Pos	Band	EA/OA	Gene
rs740681	19	3,533,276	19p13.3	G/A	FZR1	4,220	819.3	1.66 × 10–6	0.00435	0.000224	0.05	1.24 × 10–5
rs2272303	12	53,105,721	12q13.13	A/G	SOAT2	−5514	1,412	3.20 × 10–4	−0.0096	0.0027	3.10 × 10–4	1.46 × 10–5
rs11170413	12	53,105,750	12q13.13	A/G	SOAT2	−5476	1,430	3.99 × 10–4	−0.0097	0.0027	2.90 × 10–4	1.67 × 10–5
rs2272302	12	53,105,703	12q13.13	C/G	SOAT2	−4933	1,387	9.13 × 10–4	−0.0096	0.0027	3.10 × 10–4	3.67 × 10–5

FZR1: X58380.1 (GenBank).

SOAT2: X51352.1 (GenBank).

Abbreviations: Beta, regression coefficient; Chr, chromosome; EA/OA, effect allele/non‐effect allele; Pos, position; SE, standard error of Beta.

rs740681 is located in the FZR1 (discovery p = 1.66 × 10–6, replication p = .05), rs2272303, rs11170413, and rs2272302 are located in the SOAT2 (discovery p = 3.99 × 10–4, 2.90 × 10–4, and 9.13 × 10–4, respectively; replication p = 3.10 × 10–4, 2.90 × 10–4, and 3.10 × 10–4, respectively).

Meta‐analysis p‐values are ranging from 1.24 × 10–5 to 3.67 × 10–5.

Manhattan plot of the discovery cohort is displayed in Figure 2. Regional plot of the genes FZR1 and SOAT2 are drawn by LocusZoom (Figure 3a,b, respectively).

Figure 2

Manhattan plot of the discovery sample

Figure 3

Regional plots. Regional plots of the discovery samples around genes FZR1 (a) and SOAT2 (b) are presented

We annotated each of the four replicated SNPs and their neighbor SNPs (LD r 2 ≥ 0.8) with HaploReg. rs740681 is an intron INDEL variant in FZR1. It has no neighbor variant with strong LD structure. rs2272303, rs11170413, and rs2272302 are mapping to a loci: 12q13.13. The leading SNP rs1117041 has 11 SNPs in strong LD pattern, and one SNP rs2272296 has the enhancer activity marked by H3K4me1 histone mark in skeletal muscle myoblast cells, which is also significant at the discovery stage (p = 1.83 × 10–3).

From the skeletal muscle tissue summary results downloaded from the GTEx website (V7), we found that rs740681 shows an eQTL signal with DOHH in skeletal muscle (p = 1.80 × 10–14, Figure 4a). rs1170413, rs2272303, and rs2272302 show a cis‐eQTL signal with SOAT2 in Muscle‐Skeletal (p = 1.50 × 10–5, Figure 4b‐d).

Figure 4

Represents an expression quantitative trait locus (eQTL) for seven significant SNPs expression in Muscle‐Skeletal from the GTEx eQTL Browser. The Y‐axis is the rank‐normalized gene expression. rs740681 shows an eQTL signal with DOHH in Muscle‐Skeletal (p = 1.80 × 10–14, (a)). rs1170413, rs2272303, and rs2272302 show a cis‐eQTL signal with SOAT2 in Muscle‐Skeletal (p = 1.50 × 10–5, (b) (c) (d))

In the study of peripheral blood by Westra et al, rs2272302 shows a cis‐eQTL signal with gene ITGB7 (p = 1.22 × 10–17).

Gene–gene interaction network analysis of FZR1 connected them to muscle‐related genes such as MYOD1, MYOG, and so on (Figure 5).

Figure 5

Protein–protein interactions (PPI) network of FZR1 connected them to muscle‐related genes such as MYOD1, MYOG, and so on. Proteins in the interaction network were represented with nodes, while the interaction between any two proteins therein was represented with an edge. Line color indicates the type of interaction evidence including known interactions, predicted interactions, and other. These interactions contain direct (physical) and indirect (functional) interactions, derived from numerous sources such as experimental repositories and computational prediction methods. The figure was plotted by STRING

Two genes have DXA scan results in the IMPC database. Heterozygous FZR1‐deficiency mice have increased bone mineral content (p = .01).

DISCUSSION

In this study, we have performed WES study in 101 unrelated Chinese samples and replicated in two GWAS analyses. Combining the evidence from both the discovery and the replication samples, we have identified novel sarcopenia risk genes.

Previous studies implicate the ubiquitin‐proteasome system (UPS) in myogenic differentiation through regulating cell cycle progression and modulating myogenic factors such as MYOD and MYF5 (Li, Wu, & Wan, 2007). Certain ubiquitin protein ligases regulate muscle differentiation, such as SKP1‐CUL1‐F‐box protein (SCF) and anaphase‐promoting complex (APC), which have been suggested to govern terminal muscle differentiation. SCF have been suggested to govern MYOD protein degradation, while APC has been suggested to control MYF5 proteolysis during the lineage of muscle differentiation (Gardrat, Montel, Raymond, & Azanza, 1997; Lindon, Albagli, Domeyne, Montarras, & Pinset, 2000; Tintignac et al., 2005).

FZR1 (fizzy and cell division cycle 20 related 1) is a substrate‐specific adapter for the anaphase promotion complex/cyclosome (APC/C) E3 ubiquitin‐protein ligase complex and activation of APC/C in anaphase and late stage. The APC/C remains active in degrading substrates to ensure that positive regulators of the cell cycle do not accumulate prematurely. FZR1 and CDH1 (cadherin 1) share the CDH1 symbol. CDH1 is sometimes used as an abbreviation for “CDC20 homolog 1”, which is an alternate name for the FZR1. CDH1 is expressed in skeletal muscle cells (Gieffers, Peters, Kramer, Dotti, & Peters, 1999; Wan & Kirschner, 2001; Yoshiyuki, Judith, Takahiko, Shirin, & Azad, 2004). Results suggest that multiple CDH1 may temporally and spatially regulate APC activity both within and outside of the cell cycle. CDH1 may affect muscle cells (Wenqi, George, & Yong, 2007). CDH1‐APC is an important ubiquitin E3 ligase that regulates muscle differentiation by coordinating cell cycle progression and initiating a myogenic differentiation program.

SOAT2 (sterol O‐acyltransferase 2) is a member of a small family of acyl coenzyme A: cholesterol acyltransferases. The gene encodes a membrane‐bound enzyme localized in the endoplasmic reticulum that produces intracellular cholesterol esters from long‐chain fatty acyl COA and cholesterol. The cholesterol esters are then stored as cytoplasmic lipid droplets inside the cell. Several alternatively spliced transcript variants of this gene have been described, but their full‐length nature is not known. It biased expression in small intestine (RPKM 4.5), duodenum. The gene has been found to be a major regulator of cholesterol metabolism and absorption factor Buhman et al., 2000). Impaired absorption of cholesterol is associated with high body mass index (BMI) and obesity (Miettinen & Gylling, 2000; Simonen, Gylling, Howard, & Miettinen, 2000).

Conclusion: By performing WES, we identified novel sarcopenia risk variants that may play roles in skeletal muscle metabolism. Our findings may provide useful information for understanding the mechanisms of sarcopenia.

CONFLICT OF INTEREST

The authors declare that they have no conflict of interest.

AUTHOR CONTRIBUTIONS

Conceived and designed the experiments: H.W.D. and J.Y.W. Performed the experiments: S.R., L.Z., and X.H. Analyzed the data: L.Z., X.H, Z.X.J., Y.X.Z., Y.L., G.S.G Q.T., and Y.H.Z. Literature search: X.H., Y.X.Z., B.L.L, and Y.L., Wrote the paper: S.R., X.H., J.Y.W., and H.W.D. All authors reviewed and approved the manuscript.

ETHICS APPROVAL AND CONSENT TO PARTICIPATE

The study was approved by the Institutional Review Boards of Creighton University, University of Missouri‐Kansas City. Before entering the study, signed informed consent documents were obtained from all the participants.

CONSENT FOR PUBLICATION

No applicable.