Transcriptome-wide association study identifies novel genes associated with child body mineral density and lean body mass

doi:10.21203/rs.3.rs-1630560/v1

Download PDF

Research Article

Transcriptome-wide association study identifies novel genes associated with child body mineral density and lean body mass

https://doi.org/10.21203/rs.3.rs-1630560/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 27 Dec, 2022

Read the published version in Endocrine →

You are reading this latest preprint version

Objective

To identify novel candidate genes, the expression of which are associated with child body mineral density (BMD) and body lean mass (LM).

Methods

The tissue specific transcriptome-wide association (TWAS) was conducted utilizing a large scale GWAS dataset of BMD and LM, totally involving 10,414 participants. The measurement of BMD and LM phenotypes was made by TB-DXA scans. TWAS was conducted by FUSION software. The reference panels of muscle skeleton (MS), peripheral blood (NBL) and whole blood (YBL) were used for TWAS analysis. Functional enrichment and protein-protein interaction analyses of the genes identified by TWAS were performed by the online tool of Metascape (http://metascape.org).

Results

For BMD, we detected 120, 86 and 174 genes with TWAS P value < 0.05 for MS, NBL and YBL respectively. We also identified 4 common BMD associated genes shared by MS, NBL and YBL, including ZSWIM7, FAM118A, CRIPAK and ZNF641. For LM, we detected 145, 94 and 208 genes with TWAS P value < 0.05 for MS, NBL and YBL respectively. And 3 genes were detected in all of MS, NBL and YBL, including IGHMBP2, TRIT1 and LTA4H. GO enrichment analysis of BMD associated genes detected 200 GO terms, such as steroid hormone mediated signaling pathway (LogP = -3.13), autophagy (LogP = -2.65), cellular response to glucocorticoid stimulus (LogP = -2.95). GO enrichment analysis of LM associated genes detected 287 GO terms, such as Diseases of carbohydrate metabolism (LogP = -2.43), regulation of lipid localization (LogP = -2.80) and lipid storage (LogP = -3.55).

Conclusion

This study identified several candidate genes for child BMD and LM, providing novel clues for revealing the genetic mechanism underlying the development of child BMD and LM.

transcriptome-wide association study (TWAS)

genome wide association study (GWAS)

body mineral density

body lean mass

expression-trait association

Body mineral density (BMD) is the amount of bone mineral in bone tissue, often measured by the dual-energy X-ray absorptiometry (DXA). It is very widely used in clinical practice to assess the health status of bone in young people, indicating the osteoporosis and evaluating the risk of fracture. Lean mass (LM) means a composition of our body, consisting of muscle, organs, and bone theoretically, and can also be obtained by DXA. It is often thought to be a good proxy of skeletal muscle. Both of bone and skeletal muscle are the important components of musculoskeletal system, associating with the locomotion of our body.

BMD is a complex trait and its heritability was estimated to 46%-84%[1]. In the past decades, several genome-wide association studies have been conducted to identify the susceptibility genes of BMD [2–10]. Over 60 different loci were identified to be robustly associated with BMD at different skeletal sites, including PTCH1, EN1, WNT16, ESR1, LGR4 and so on[11]. Meanwhile, the heritability of LM was estimated to around 52%[11]. However, there were limited genome-wide association studies of LM compared with BMD [12–14]. The associated genes included PRDM16, GLYAT et al. Moreover, there existing a genetic correlation between BMD and LM ranging from 69%-88% depending on the skeletal sites[15]. It was thought that enhancing the health status of bone, especially maximizing the peak bone mass in early life would help reduce the risk of osteoporosis later in adulthood[16]. So to understand the regulation mechanism of child bone phenotypes is critically important. But few efforts have been paid into the study.

In the classic GWASs, vast association signals have been detected between genetic variants and traits or phenotypes. In most cases, the association does not mean the real causal relationship partly for the mechanistic steps between genetic variants and traits or phenotypes cannot be taken into consideration. Furthermore, several evidences show that a substantial proportion of risk variants exerts their influence on traits by modulating the expression levels of the target genes (for example, in the case of eQTLs)[17, 18]. So in such context, the notion of expression-trait association is raised. But it has been hindered for the difficulties in the specimen collection and the high cost of genotype, phenotype and gene expression measurement meanwhile. For addressing these issues, the transcriptome-wide association study (TWAS) is proposed, which can identify the genes whose expression is associated with the complex traits based on the GWAS summary data[19]. Without directly measuring the gene expression, this approach leveraged a relative small set of reference panel with both genotype and expression data to impute the gene expression level in a large-scale GWAS using the genotype data. And then, the imputed gene expression was used to identify the expression-trait association. Through this approach, several expression-trait associations were found in various traits or diseases [19–23]. For example, Gusev, A et al. detected 157 TWAS-significant genes in schizophrenia GWAS and finally found one of the genes, mapk3, showed a significant effect on neurodevelopmental phenotypes in zebrafish[20].

For obtaining more interpretable biological unit of BMD and LM, we conducted this tissue specific TWAS based on a large GWAS dataset of BMD and LM. The reference panels in different tissues/cells including muscle skeleton (MS), peripheral blood (NBL) and whole blood (YBL) were used for the imputation. Several significant expression-trait association signals were detected in both BMD and LM. Moreover, an enrichment analysis was performed to explore the functional annotation.

GWAS summary data

The GWAS data used in our analysis was extracted from a previously published study[11]. This study enrolled 10,414 participants from four pediatric cohorts, including the Generation R Study, the Avon Longitudinal Study of Parents and their Children (ALSPAC), the Bone Mineral Density in Childhood Study (BMDCS), and the Copenhagen Prospective Studies on Asthma in Childhood (COPSAC) cohort. Average ages of participants from four cohorts varied in a range of 6.21–9.94. The measurement of BMD and LM phenotypes were made by TB-DXA scans. DNA samples were genotyped using the Illumina platform. The SNPs with MAF < 0.05 were excluded. Imputation to the CEU panel of the HapMap Phase II (build 36 release 22) reference was performed. BMD and LM measurements were adjusted by age, gender, height, fat percent (TB-FM/weight), and study specific covariates (genetic principal components and measurement center). And the GCTA software package as well as LD-score regression methodology was used to assess the SNP heritability.

The transcriptome-wide association study (TWAS)

TWAS was performed using the FUSION software based on the pediatric GWAS summary data (http://gusevlab.org/projects/fusion/). The software could provide an approach to identify the association between the complex traits and intermediate phenotype (gene expression) without directly measuring the expression levels[19]. Reference panels of different tissues was used to impute the cis genetic component expression of SNP genotype data. In our study, the gene expression weights of muscle skeleton (MS), peripheral blood (NBL) and whole blood (YBL), which was driven from the FUSION, were used for TWAS analysis. A P value was calculated by TWAS for each analyzed gene.

Functional annotation and enrichment analysis

For better interpreting the associated genes found in the TWAS, the online tool, Metascape, was used to make further functional enrichment analysis and protein-protein interaction analysis of the genes identified by TWAS (http://metascape.org). It could apply a standard accumulative hypergeometric statistical test to identify ontology terms based on the resources of Canonical Pathway (MSigDB), Hallmark Gene Sets (MSigDB), KEGG Pathway and GO Bilogical Processes. The construction of PPI-network (protein-protein interaction network) and associated module analysis were based on GO enrichment analysis using the plugin Molecular Complex Detection (MCODE). MCODE algorithm was then applied to this network to identify neighborhoods where proteins are densely connected.

The associated genes found in the TWAS

For different tissue/cell of MS, NBL and YBL, 2976, 2455 and 4701 genes were finally included in the expression-train association analysis (seen in Table 1). For child BMD phenotype, TWAS identified 120, 86 and 174 genes with TWAS P value < 0.05 for MS, NBL and YBL, respectively. For child LM, we detected 145, 94 and 208 genes with TWAS P value < 0.05f or MS, NBL and YBL, respectively.

Table 1

**The number of associated with BMD and LM identified by the TWAS**
Tissue/cell	N^#	n^* for BMD	n^* for LM
MS	2976	120	145
NBL	2455	86	94
YBL	4701	174	208
Note: N means the number of genes finally included in the analysis of TWAS.
n* means the number of significantly associated genes identified in the TWAS.

Moreover, a venn diagram was made to found the overlapping significantly associated genes in three different tissues/cells for BMD and LM (See in Fig. 1). For BMD, we found 4 common genes shared by MS, NBL and YBL, including ZSWIM7 (P_MS= 1.85×10^− 2, P_NBL=3.51×10^− 2, P_YBL=2.27×10^− 2), FAM118A(P_MS=7.87×10^− 4, P_NBL=6.14×10^− 3, P_YBL=6.02×10^− 3), 3CRIPAK(P_MS=4.58×10^− 2, P_NBL=1.04×10^− 2, P_YBL=2.62×10^− 4) and ZNF641(P_MS=3.66×10^− 2, P_NBL=1.19×10^− 3, P_YBL=1.23×10^− 3). For LM, we found 3 common genes detected in all of MS, NBL and YBL, including IGHMBP2 (P_MS= 2.12×10^− 2, P_NBL=2.39×10^− 2, P_YBL=3.34×10^− 2), TRIT1(P_MS= 2.21×10^− 2, P_NBL= 2.50×10^− 2, P_YBL=3.39×10^− 2) and LTA4H(P_MS=1.91×10^− 2, P_NBL=1.72×10^− 2, P_YBL=1.78×10^− 2). Table 2 presented the detailed information of the 7 common genes.

Table 2

**The detailed information of all the seven common associated genes**
Gene Symbol	HSQ			BEST.GWAS.ID	NSNP	TWAS P value
Gene Symbol	MS	NBL	YBL	BEST.GWAS.ID	NSNP	MS	NBL	YBL
BMD
ZSWIM7	0.22	0.02	0.63	rs12945204	272	1.85×10^− 2	3.51×10^− 2	2.27×10^− 2
FAM118A	0.32	0.41	0.49	rs3827393	617	7.87×10^− 4	6.14×10^− 3	6.02×10^− 3
CRIPAK	0.25	0.13	0.35	rs2335937	326	4.58×10^− 2	1.04×10^− 2	2.62×10^− 4
ZNF641	0.19	0.27	0.43	rs12306165	473	3.66×10^− 2	1.19×10^− 3	1.23×10^− 3
LM
IGHMBP2	0.23	0.02	0.23	rs11228269	406	2.12×10^− 2	2.39×10^− 2	3.34×10^− 2
TRIT1	0.16	0.03	0.26	rs7415236	379	2.21×10^− 2	2.50×10^− 2	3.39×10^− 2
LTA4H	0.20	0.11	0.56	rs1108405	550	1.91×10^− 2	1.72×10^− 2	1.78×10^− 2
Note: HSQ means heritability of genes, BEST.GWAS.ID means rsID of the most significant GWAS SNP in locus, NSNP means number of SNPs in the locus.

Functional enrichment and PPI analysis results

The heatmap of the enriched terms was shown in Fig. 2. For BMD, GO enrichment analysis detected 200 GO terms, such as steroid hormone mediated signaling pathway (GO: 0043401, LogP=-3.13), autophagy (GO: 0006914, LogP=-2.65), regulation of cytokine secretion (GO: 0050707, LogP=-3.43) and cellular response to glucocorticoid stimulus (GO: 0071385, LogP=-2.95). For LM, GO enrichment analysis detected 287 GO terms, such as Diseases of carbohydrate metabolism (R-HSA-5663084, LogP=-2.43), regulation of lipid localization (GO: 1905952, LogP=-2.80), lipid storage (GO: 0019915, LogP=-3.55) and mitochondrial RNA processing (GO: 0000963, LogP=-3.41). Protein-protein interaction enrichment analysis was also carried out with the following databases: BioGrid7, InWeb_IM8, OmniPath9. The MCODE networks identified for individual gene lists were pooled and shown in Fig. 3.

BMD and LM are the complex traits which can be influenced by multiple genetic factors. Although various genes have been identified to be associated with the two traits respectively by previous studies, limited mechanistic clues can be obtained from these studies. It is known that the genetic variants can exert influence on the traits by regulating the gene expression level. So we conducted this tissue specific transcriptome-wide association study to investigate the expression-trait association for BMD and LM. The databases with both expression and genotype measurement in the tissues/cells including muscle skeleton, peripheral blood and whole blood were used as reference panel for the expression imputation.

IGHMBP2 (immunoglobulin mu DNA binding protein 2) is one of the significant associated genes with LM across LM, YBL and NBL. It locates in 11q13.2 and can encode a kind of helicase to switch immunoglobulin by binding to a specific DNA sequence. Mutation of this gene is thought to be involved in the spinal muscle atrophy with respiratory distress type 1(SMARD1), which can cause muscle weakness and respiratory failure typically beginning in infancy[24]. The deficiency of IGHMBP2 protein can lead to the degeneration of muscle cell nuclei [25]. But how this gene affects the phenotype of LM remains unclear and need more studies.

MTHFR is another notable gene associated with LM for MS and YBL. The protein encoded by this gene is a kind of enzyme in the methyl cycle which can catalyze the conversion of 5, 10-methylenetetrahydrofolate to 5-methyltetrahydrofolate. In a previous study conducted by Liu[26], the MTHFR gene polymorphism was found to be associated with body lean mass but not fat body mass, which was consistent with our study results.

CRIPAK was identified to be significantly associated with BMD in all three tissues. It had a role in the modulation of Pak1-mediated estrogen receptor transactivation by negatively regulating Pak1[27]. The estrogen stimulation of cells could enhance the CRIPAK expression and promote its co-localization with estrogen receptor in the nuclear compartment [27]. Estrogen plays an important role in the growth and maturation of bone. At cellular level in bone estrogen can decrease the cell number and reduce the amount of active remodeling units by inhibiting differentiation of osteoclasts[28]. So CRIPAK might exert an influence on the BMD by the estrogen related pathway.

MSH3 was also found to be significantly associated with BMD in MS and YBL. The protein encoded by this gene can bind to MSH2 and form MutS beta, which belongs to the post-replicative DNA mismatch repair system. Previous studies identified that some SNP loci on this gene (such as rs2035256, rs33013) was significantly associated with BMD and have strong cis-effects on gene expression in human primary osteoblast [29, 30]. But for the lack of further studies of its function, how it influences the BMD phenotype remains unknown.

The enrichment analysis results showed that the significant associated genes with BMD were enriched in the hormone-related categories, such as steroid hormone mediated signaling pathway (GO: 043401), regulation of growth hormone secretion (GO: 0060123) and cellular response to glucocorticoid stimulus (GO: 0071385). The steroid hormone is a kind of hormone making of steroid compound, mainly secreted by tissues including adrenal cortex, testes and ovaries. It can be divided into five types according to their receptors, such as glucocorticoids, mineralocorticoids, androgens, estrogens, and progestogens. Estrogens are known to be the powerful regulator of bone metabolism. The loss of estrogen from the ovarian after menopause is always accompanied with the decline of bone mineral density, which lead to higher osteoporosis rate in women than men. The glucocorticoid is associated with a decreased bone mineral density and impaired bone microarchitecture parameters [31]. Growth hormone is a peptide hormone which can stimulate cell reproduction and cell regeneration in humans and play an important role in human development. These findings seem consistent with the previous knowledge of BMD.

For LM, it seems that the associated genes are more likely to be enriched in the metabolism categories of several substances, including the heterocycle catabolic metabolism (GO: 0046700), carbohydrate metabolism (R-HAS-5663084), regulation of lipid localization (GO: 1905952), tetrapyrrole metabolism (GO: 0033013), purine metabolism (hsa00230), lipid storage (GO: 0019915) and monocarboxylic acid metabolism process (GO: 0032787). These are all important nutrition substances of our body. Carbohydrate and lipid are two kinds of essential macronutrients for the organism. A heterocycle is a cyclic compound containing atoms of at least two different elements in its ring, such as purine. Purine is the basic components of nucleic acid. Tetrapyrroles are a group of chemical compounds containing four pyrrole or pyrrole-like rings and the core parts of some compounds with crucial biochemical roles in the living system, such as hemoglobin. And monocarboxylic acids are the essential resources for synthesizing amino acid. Compared with BMD, LM is more likely to be regulated by the biochemical metabolism process.

The TWAS analysis approach developed by Alexander G et al. was adopted in this study[19]. It can identify genes whose expression is significantly associated with complex traits in individuals without directly measuring the expression levels based on the GWAS summary data[19]. There are several potential advantages by using this approach. First of all, the expression-trait association can provide more interpretable clues for the genetic study while the GWAS often obtained associated locus lying in the LD with multiple significant SNPs which may not in the genes. Moreover, unlike the analysis focusing on eQTL and SNP association, TWAS can combine full cis-SNP signals, no matter they are significant or not, to make an expression imputation. Finally, it can also avoid confounding from environment differences caused by the traits.

There are also some issues to be addressed about the limitations of this study. First of all, as the gene expression is imputed based on the reference panels, there is a possibility that the results can be influenced by the quality and sample size of the reference data. A larger sample size and more available tissues’ datasets can mitigate such impact. Second, the number of imputed genes also depends on the training data. So it is limited in some degree. Finally, only the expression mediated regulation modes are considered by this analysis approach although there are lots of other ways by which the SNPs influence the traits.

This study conducted a tissue specific TWAS for the phenotypes of BMD and LM based on the previous GWAS datasets. The expressions of several genes were identified to be associated with BMD and LM respectively. By performing gene enrichment analysis, we found a different regulation mechanism of the two traits. This study could provide novel clues for the study of BMD and LM and outline a systematic approach to identify functional mediators of complex disease.

Ethics approval and consent to participate

The source of the data was a publicly available data base and no human participants were involved, hence ethical parameters are not applicable.

Consent for publication

Not applicable

Availability of data and supporting materials

The GWAS summary dataset was extracted from a previously published study (Bivariate genome-wide association analysis implicates pleiotropic effects at the SREBF1/TOM1L2 locus on bone mineral density and lean mass in children)

Competing interests

The authors declare that they have no competing interests.

Funding

This work was supported by the National Natural Science Foundation of China (grant number 81974347 and 81802210); the Department of Science and Technology of Sichuan Province (grant number 2021YFS0122). Financial support had no impact on the outcomes of this study.

Acknowledgements

Not applicable

Author details

Orthopedic Research Institute, Department of Orthopedics, Sichuan University West China Hospital, 37# Guoxue Road, Chengdu, 610041, Sichuan Province, People’s Republic of China

N.K. Arden et al., The heritability of bone mineral density, ultrasound of the calcaneus and hip axis length: A study of postmenopausal twins. J. Bone Mineral Res. 11(4), 530–534 (2010)
U. Styrkarsdottir et al., Sequence variants in the PTCH1 gene associate with spine bone mineral density and osteoporotic fractures. Nat. Commun. 7, 10129 (2016)
H.F. Zheng et al., Whole-genome sequencing identifies EN1 as a determinant of bone density and fracture. Nature 526(7571), 112–117 (2015)
R. Fernando et al., Twenty bone mineral density loci identified by large-scale meta-analysis of genome-wide association studies. Nat. Genet. 41(11), 1199–1206 (2009)
J.B. Richards et al., Bone mineral density, osteoporosis, and osteoporotic fractures: a genome-wide association study. Lancet 371(9623), 1505–1512 (2008)
K. Estrada et al., Genome-wide meta-analysis identifies 56 bone mineral density loci and reveals 14 loci associated with risk of fracture. Nat. Genet. 44(5), 491 (2012)
D.L. Koller et al., META-ANALYSIS OF GENOME-WIDE STUDIES IDENTIFIES WNT16 AND ESR1 SNPS ASSOCIATED WITH BONE MINERAL DENSITY IN PREMENOPAUSAL WOMEN. J. Bone Mineral Res. 28(3), 547–558 (2013)
U. Styrkarsdottir et al., Nonsense mutation in the LGR4 gene is associated with several human diseases and other traits. Nature 497(7450), 517–520 (2013)
U. Styrkarsdottir et al., Two Rare Mutations in the COL1A2 Gene Associate With Low Bone Mineral Density and Fractures in Iceland. J. Bone Mineral Res. 31(1), 173–179 (2016)
L. Zhang et al., Multistage genome-wide association meta-analyses identified two new loci for bone mineral density. Hum. Mol. Genet. 23(7), 1923–1933 (2014)
M.G.E. Al, Bivariate genome-wide association analysis implicates pleiotropic effects at the SREBF1/TOM1L2 locus on bone mineral density and lean mass in children. Nature Communications, 2017. -(-)
R. Hai et al., Genome-wide association study of copy number variation identified gremlin1 as a candidate gene for lean body mass. J. Hum. Genet. 57(1), 33–37 (2012)
T. Urano et al., Large-scale analysis reveals a functional single‐nucleotide polymorphism in the 5′‐flanking region of PRDM16 gene associated with lean body mass. Aging Cell,13,4(2014-05-23), 2014. 13(4): p. 739–743
Y.F. Guo et al., Suggestion of GLYAT gene underlying variation of bone size and body lean mass as revealed by a bivariate genome-wide association study. Hum. Genet. 132(2), 189–199 (2013)
L.H. Bogl et al., An investigation into the relationship between soft tissue body composition and bone mineral density in a young adult twin sample. J. Bone Mineral Res. Official J. Am. Soc. Bone Mineral Res. 26(1), 79–87 (2011)
N.H. Golden, S.A. Abrams, Optimizing bone health in children and adolescents. Pediatrics 134(4), e1229–e1243 (2014)
N. Mancuso et al., Integrating Gene Expression with Summary Association Statistics to Identify Genes Associated with 30 Complex Traits. Am. J. Hum. Genet. 100(3), 473 (2017)
D.L. Nicolae et al., Trait-associated SNPs are more likely to be eQTLs: annotation to enhance discovery from GWAS. PLoS Genet. 6(4), e1000888 (2010)
A. Gusev, A. Ko, Integrative approaches for large-scale transcriptome-wide association studies. 2016. 48(3): p. 245–52
A. Gusev et al., Transcriptome-wide association study of schizophrenia and chromatin activity yields mechanistic disease insights. Nat. Genet. 50(4), 538–548 (2018)
S. Thériault et al., A transcriptome-wide association study identifies PALMD as a susceptibility gene for calcific aortic valve stenosis. Nature Communications, 2018. 9(1)
T. Jiang et al., Transcriptome-wide association study revealed two novel genes associated with nonobstructive azoospermia in a Chinese population. Fertility & Sterility 108(6), 1056–1062 (2017)
H. Lin et al., Transcriptome-wide association study of inflammatory biologic age. Aging 9(11), 2288–2301 (2017)
P. Francesca et al., The wide spectrum of clinical phenotypes of spinal muscular atrophy with respiratory distress type 1: a systematic review. J. Neurol. Sci. 346(1–2), 35–42 (2014)
M. Jedrzejowska et al., Severe phenotypes of SMARD1 associated with novel mutations of the IGHMBP2 gene and nuclear degeneration of muscle and Schwann cells. Eur. J. Paediatr. Neurol. 18(2), 183–192 (2014)
X. Liu et al., The MTHFR gene polymorphism is associated with lean body mass but not fat body mass. Hum. Genet. 123(2), 189–196 (2008)
A.H. Talukder, Q. Meng, R. Kumar, CRIPak, a novel endogenous Pak1 inhibitor. Oncogene 25(9), 1311–1319 (2006)
H.K. Vaananen, P.L. Harkonen, Estrogen and bone metabolism. Maturitas 23 Suppl, S65–S69 (1996)
U. Styrkarsdottir et al., Multiple genetic loci for bone mineral density and fractures. N. Engl. J. Med. 358(22), 2355 (2008)
E. Grundberg et al., Population genomics in a disease targeted primary cell model. Genome Res. 19(11), 1942–1952 (2009)
S. Adhikary et al., Dietary flavonoid kaempferol inhibits glucocorticoid-induced bone loss by promoting osteoblast survival. Nutrition (Burbank, Los Angeles County, Calif.), 2018. 53: p. 64–76

Download PDF

Journal Publication

published 27 Dec, 2022

Read the published version in Endocrine →

Reviewers agreed at journal
01 Jun, 2022
Reviews received at journal
09 May, 2022
Reviewers invited by journal
07 May, 2022
Editor assigned by journal
06 May, 2022
First submitted to journal
06 May, 2022

You are reading this latest preprint version

Transcriptome-wide association study identifies novel genes associated with child body mineral density and lean body mass

Status:

Journal Publication

Version 1

Abstract

Objective

Methods

Results

Conclusion

Figures

Introduction

Materials And Methods

GWAS summary data

The transcriptome-wide association study (TWAS)

Functional annotation and enrichment analysis

Results

The associated genes found in the TWAS

Functional enrichment and PPI analysis results

Discussion

Conclusion

Declarations

References

Status:

Journal Publication

Version 1