Martínez-Gil N, Patiño J, Ugartondo N, Grinberg D, Balcells S

Department of Genetics, Microbiology and Statistics. Biology Faculty. University of Barcelona, CIBERER, IBUB, IRSJD. Barcelona (Spain)

DOI: 10.4321/S1889-836X2021000300003

Study awarded a FEIOMM research grant 2019.

Summary
Introduction: WNT16 is an important gene in bone homeostasis, found in a very complex locus that also includes neighboring genes: ING3, FAM3C, and CPED1. In addition to the clear role of WNT16 in determining bone mineral density (BMD), evidence has also been found for the importance of these three neighboring genes in bone metabolism. Therefore, it remains to be clarified whether the variants in WNT16 associated with BMD carry out their effect on WNT16 or if they do so by modifying the expression of these neighboring genes.
Material and methods: We have determined the expression levels of CPED1 and FAM3C in primary osteoblasts and we have verified whether WNT16 variants behave as loci of quantitative expression traits (expression quantitative trait loci; eQTL) of these genes.
Results: The amino acid change variant rs2908004 in WNT16 acts as the eQTL of FAM3C in primary osteoblasts under the dominant model hypothesis.
Discussion: It is possible that the effect of this variant on BMD is due to the modification of the expression levels of FAM3C in addition to or instead of a direct effect of the mutant WNT16 protein resulting from the amino acid change.

Key words: Bone mineral density, WNT16, osteoporosis, transcription.

Introduction
WNT16 is a ligand of the Wnt pathway that has been extensively studied for its importance in regulating bone homeostasis. This has been confirmed with the phenotype of knock-out (KO) mice and conditional KO mice in osteoblasts (cKO), which show spontaneous fractures due to low cortical bone mineral density (BMD), low bone strength and high cortical porosity, keeping the volume of trabecular bone unchanged [1-4]. On the contrary, Wnt16 overexpression in osteoblasts and osteocytes produces an increase in BMD and bone strength in both trabecular and cortical bone [5-7]. Despite this, the precise mechanism by which WNT16 acts is not known and different studies indicate that the effect on canonical and non-canonical Wnt pathways could be tissue specific [1,8-11]. In bone, WNT16 is expressed mainly by osteoblasts and carries out its function both by stimulating bone formation and by inhibiting its resorption indirectly through osteoprotegerin (OPG) or directly affecting the osteoclasts differentiation [1,12].
Several genome-wide association studies (GWAS) have shown an association between the locus containing WNT16 and various skeletal phenotypes, including BMD and risk of fractures [2,3,13-26]. WNT16 is found in a very complex locus, where several genes in the region show an important role in bone metabolism. The genes ING3 and CPED1 at 5 ‘and FAM3C at 3’ of WNT16 belong to this locus (figure 1).
ING3 (Inhibitor Of Growth Family Member 3) is responsible for the regulation of chromatin, since it is part of the NuA4 histone acetyltransferase (HAT) complex that recognizes the trimethylated form of lysine 4 of histone H3 [27]. Other functions unrelated to chromatin regulation include apoptosis promotion, DNA repair, and modulation of cell mobility [27]. ING3 is expressed in a multitude of tissues, especially in those with a higher proportion of cell growth, bone being one of those with the highest expression of ING3 [27]. The in vitro cellular model of ING3 KO mesenchymal cells shows involvement of osteoblastogenesis and stimulation of adipogenic differentiation [28].
No specific function is known for CPED1 (Cadherin Like And PC-Esterase Domain Containing 1) in humans or mice. In the latter, it is uniformly expressed in a variety of solid tissues, including bone, although it is not detected in the RAW264.7 cell line or in circulating leukocytes [29]. Furthermore, CPED1 presents different isoforms due to alternative splicing and three active promoter regions during osteogenic differentiation, indicating a complex regulation during differentiation [29].
FAM3C (Family of sequence similarity 3c) is a cytokine-like growth factor expressed in a multitude of tissues [30], which plays a very important role in the epithelial-mesenchymal transition (EMT) and subsequent metastasis during cancer progression [31]. Its relationship with bone metabolism has been confirmed with the KO mouse model, which presents alterations in the cortical and trabecular structure, with an increase in cortical BMD, which generates a decrease in bone resistance [30]. In in vitro studies, it was found that the mesenchymal cells extracted from these KO mice showed accelerated osteoblastogenesis [31].
The relationship of the genes present at the ‘ING3-CPED1-WNT16-FAM3C’ locus with bone metabolism and their repeated association with BMD raises the question of whether there is a single causative gene and, if so, which is it, or if instead, all genes are contributing to the phenotype. To this end, in the present work we have determined whether those WNT16 variants associated with BMD in a previous work by our group [32] are found to modify the expression of neighboring genes CPED1 and FAM3C.

Material and methods
Cell culture
Human primary osteoblasts (hOB) were used for the loci assays that determine quantitative differences in gene expression (expression quantitative trait loci; eQTL). The hOBs were obtained from trabecular bone fragments discarded from knee replacement operations performed on women with osteoarthritis and who did not have any other disorder that could affect bone quality. The study was approved by the Clinical Research Ethics Committee of the Parc de Salut MAR (registration numbers: 2010/3882/I and 2013/5266/I) and was carried out in accordance with the Declaration of Helsinki, obtaining informed consent by written by all participants. The primary osteoblast culture protocol is described in Roca-Ayats et al. [33]. Briefly, bone samples were cut into small pieces and washed with phosphate buffered saline (PBS; Gibco, Life Technologies). These pieces were cultured in 140 mm plates with DMEM supplemented with 10% FBS, 1% w/s, 0.4% fungizone (Gibco, Life Technologies) and 100 µg/ml of ascorbic acid (Sigma-Aldrich). When the cells reached confluence, they were divided into three 75 cm2 flasks, one for DNA extraction, one for RNA extraction, and the third for freezing and storage. Cells at passage 2 or lower were used for all extractions.

EQTL assay
DNA was extracted from cultured hOBs using the Wizard® Genomic DNA Purification Kit (Promega), according to the manufacturer’s instructions. The concentration of the purified DNA and its quality was analyzed in a spectrophotometer (Nanodrop). The genotypes of the rs2908004, rs2707466, rs55710688 and rs142005327 variants were evaluated by Sanger sequencing using BigDye® Terminator v3.1 (Applied Biosystems) in the Genomics facilities of the CCiT of the University of Barcelona. The primers (Invitrogen, Thermo Fisher) were designed using Primer3 Input 0.4.0 (table 1). The total RNA of the cultured hOBs was extracted using the High Pure RNA Isolation kit (Roche), according to the manufacturer’s instructions and the quantification and quality of the RNA were checked using a Nanodrop spectrophotometer. RNA was reverse transcribed to cDNA using the High Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Thermo Fisher), according to the manufacturer’s specifications. RT-qPCR was carried out using UPL probes (Roche) on a LightCycler 480 Instrument II (Roche). HMBS gene expression was used as a normalization control and the relative quantification (fold change) was calculated using the second derivative method. The number and sequence of the probe used, as well as the primers used for the amplification of the CPED1, FAM3C and HMBS genes are shown in table 1.

Statistic analysis
For the statistical analysis of the eQTL, the WGassociation34 function was used in RStudio. This function performs an association analysis between a given SNP and a dependent variable (in this case the expression levels of CPED1 and FAM3C) in five different genetic inheritance models: codominant [homozygous for major allele vs. heterozygous vs. homozygous for minor allele], dominant [homozygous for major allele versus (heterozygous + homozygous for minor allele)], recessive [homozygous for minor allele versus (heterozygous + homozygous for major allele)], over-dominant [heterozygous versus (homozygous for allele major + homozygous for minor allele)] and log-additive [each allele modifies the risk by an additive amount] .

Results
Cis-eQTL analysis
The variants rs2908004, rs2707466, rs55710688 and rs142005327 of WNT16 have been described as cis-eQTLs, according to the GTEx database in different human tissues (table 2). Unfortunately, this database does not have information on any bone tissue. This is why, using our own database of human primary osteoblasts (n=45), we have tested whether these variants act as cis-eQTL of the neighboring genes of WNT16: CPED1 and FAM3C. Only the rs2908004 variant has shown a significant association with FAM3C expression levels under the dominant hypothesis (p=0.03, table 3, figure 2). In addition, the rs2908004 and rs2707466 variants show a trend towards significance with FAM3C expression levels under the codominance hypothesis (p=0.05491) and under the dominant hypothesis (p=0.06954), respectively (table 3, figure 2). The presence of the G allele (rs2908004) and the C allele (rs2707466) are associated with an increase in the expression of FAM3C (table 3, figure 2). On the contrary, we have not found a significant association or trend between the WNT16 variants analyzed and CPED1 expression levels (table 3, figure 2).

Discussion
Over the past 20 years, many works have highlighted the importance of WNT16 in bone homeostasis. WNT16 is found at position 7q31.31 along with 3 other genes also related to bone metabolism. In a previous study by our group, we found that the variants rs142005327, rs55710588, rs2908004, and rs2707466 were associated with BMD in a cohort of postmenopausal women from the Barcelona area (BARCOS) [32]. To determine whether these variants are related to an effect on WNT16 or to an effect on the expression of neighboring genes, we have verified whether they are acting as eQTL of FAM3C and CPED1. This work has allowed us to determine that the missense variant rs2908004 is acting as eQTL of the FAM3C gene under the hypothesis of a dominant model in human primary osteoblasts.
The missense variant rs2908004 (p.Gly72Arg/ p.Gly82Arg) has been associated with different bone parameters by us and others [2,32,35-38]. This amino acid change from glycine to arginine is considered tolerated and benign by the pathogenicity predictors SIFT and PolyPhen-2, so its effect on BMD could be due to its role as eQTL and not to a change in the resulting WNT16 protein.
It should be taken into account that to obtain a more robust statistical significance having a bank with a greater number of primary osteoblasts is required. Unfortunately obtaining these samples is difficult and we have only managed to enter 45 samples into our bank.
In addition, it would be interesting to determine the expression levels of other neighboring genes that may be influencing BMD such as ING3 or directly on the expression levels of WNT16, which have not been able to be quantified due to lack of RNA sample from primary osteoblasts.

Conclusion
Through this work we have determined that the variant rs2908004 of WNT16 regulates the expression levels of the neighboring gene FAM3C under the hypothesis of a dominant model. If this association is confirmed in a larger primary osteoblast bank, it would indicate that the association of this variant with BMD could be due, at least in part, to the variation of FAM3C expression.

Conflict of interests: The authors declare no conflict of interest.

Bibliography

1. Movérare-Skrtic S, Henning P, Liu X, Nagano K, Saito H, Börjesson AE, et al. Osteoblast-derived WNT16 represses osteoclastogenesis and prevents cortical bone fragility fractures. Nat Med. 2014; 20(11): 1279-88.
2. Zheng H-F, Tobias JH, Duncan E, Evans DM, Eriksson J, Paternoster L, et al. WNT16 influences bone mineral density, cortical bone thickness, bone strength, and osteoporotic fracture risk. PLoS Genet. 2012; 8(7): e1002745.
3. Medina-Gomez C, Kemp JP, Estrada K, Eriksson J, Liu J, Reppe S, et al. Meta-analysis of genome-wide scans for total body BMD in children and adults reveals allelic heterogeneity and age-specific effects at the WNT16 locus. PLoS Genet. 2012; 8(7): e1002718.
4. Wergedal JE, Kesavan C, Brommage R, Das S, Mohan S. Role of WNT16 in the regulation of periosteal bone formation in female mice. Endocrinology. 2015; 156(3): 1023-32.
5. Alam I, Alkhouli M, Gerard-O’Riley RL, Wright WB, Acton D, Gray AK, et al. Osteoblast-Specific Overexpression of Human WNT16 Increases Both Cortical and Trabecular Bone Mass and Structure in Mice. Endocrinology. 2016; 157(2): 722-36.
6. Alam I, Reilly AM, Alkhouli M, Gerard-O’Riley RL, Kasipathi C, Oakes DK, et al. Bone Mass and Strength are Significantly Improved in Mice Overexpressing Human WNT16 in Osteocytes. Calcif Tissue Int. 2017; 100(4): 361-73.
7. Movérare-Skrtic S, Wu J, Henning P, Gustafsson KL, Sjögren K, Windahl SH, et al. The bone-sparing effects of estrogen and WNT16 are independent of each other. Proc Natl Acad Sci U S A. 2015; 112(48): 14972-7.
8. Jiang Z, Von den Hoff JW, Torensma R, Meng L, Bian Z. Wnt16 is involved in intramembranous ossification and suppresses osteoblast differentiation through the Wnt/β-catenin pathway. J Cell Physiol. 2014; 229(3): 384-92.
9. Hendrickx G, Boudin E, Verbeek M, Fransen E, Mortier G, Van Hul W. WNT16 Requires Gα Subunits as Intracellular Partners for Both Its Canonical and Non-Canonical WNT Signalling Activity in Osteoblasts. Calcified Tissue International. 2020; 106(3): 294-302.
10. Shen J, Chen X, Jia H, Meyers CA, Shrestha S, Asatrian G, et al. Effects of WNT3A and WNT16 on the Osteogenic and Adipogenic Differentiation of Perivascular Stem/Stromal Cells. Vol. 24, Tissue Engineering Part A. 2018. p. 68-80.
11. Sebastian A, Hum NR, Morfin C, Murugesh DK, Loots GG. Global gene expression analysis identifies Mef2c as a potential player in Wnt16-mediated transcriptional regulation. Vol. 675, Gene. 2018. p. 312-21.
12. Kobayashi Y, Thirukonda GJ, Nakamura Y, Koide M, Yamashita T, Uehara S, et al. Wnt16 regulates osteoclast differentiation in conjunction with Wnt5a. Biochem Biophys Res Commun. 2015; 463(4): 1278-83.
13. Moayyeri A, Hsu Y-H, Karasik D, Estrada K, Xiao S-M, Nielson C, et al. Genetic determinants of heel bone properties: genome-wide association meta-analysis and replication in the GEFOS/GENOMOS consortium. Hum Mol Genet. 2014; 23(11): 3054-68.
14. Kemp JP, M orris JA, Medina-Gomez C, Forgetta V, Warrington NM, Youlten SE, et al. Identification of 153 new loci associated with heel bone mineral density and functional involvement of GPC6 in osteoporosis. Nat Genet. 2017; 49(10): 1468-75.
15. Medina-Gomez C, Kemp JP, Trajanoska K, Luan J ’an, Chesi A, Ahluwalia TS, et al. Life-Course Genome-wide Association Study Meta-analysis of Total Body BMD and Assessment of Age-Specific Effects. Am J Hum Genet. 2018; 102(1): 88-102.
16. Estrada K, Styrkarsdottir U, Evangelou E, Hsu Y-H, Duncan EL, Ntzani EE, et al. Genome-wide meta-analysis identifies 56 bone mineral density loci and reveals 14 loci associated with risk of fracture. Nat Genet. 2012; 44(5): 491-501.
17. Trajanoska K, Schoufour JD, de Jonge EAL, Kieboom BCT, Mulder M, Stricker BH, et al. Fracture incidence and secular trends between 1989 and 2013 in a population based cohort: The Rotterdam Study. Bone. 2018; 114: 116-24.
18. Morris JA, Kemp JP, Youlten SE, Laurent L, Logan JG, Chai RC, et al. An atlas of genetic influences on osteoporosis in humans and mice. Nat Genet. 2019; 51(2): 258-66.
19. Kichaev G, Bhatia G, Loh P-R, Gazal S, Burch K, Freund MK, et al. Leveraging Polygenic Functional Enrichment to Improve GWAS Power. Am J Hum Genet. 2019; 104(1): 65-75.
20. Chesi A, Mitchell JA, Kalkwarf HJ, Bradfield JP, Lappe JM, McCormack SE, et al. A trans-ethnic genome-wide association study identifies gender-specific loci influencing pediatric aBMD and BMC at the distal radius. Hum Mol Genet. 2015; 24(17): 5053-9.
21. Koller DL, Zheng H-F, Karasik D, Yerges-Armstrong L, Liu C-T, McGuigan F, et al. Meta-analysis of genome-wide studies identifies WNT16 and ESR1 SNPs associated with bone mineral density in premenopausal women. J Bone Miner Res. 2013; 28(3): 547-58.
22. Pei Y-F, Hu W-Z, Yan M-W, Li C-W, Liu L, Yang X-L, et al. Joint study of two genome-wide association meta-analyses identified 20p12.1 and 20q13.33 for bone mineral density. Bone. 2018; 110: 378-85.
23. Zhang L, Choi HJ, Estrada K, Leo PJ, Li J, Pei Y-F, et al. Multistage genome-wide association meta-analyses identified two new loci for bone mineral density. Hum Mol Genet. 2014; 23(7): 1923-33.
24. Wang H, Zhang F, Zeng J, Wu Y, Kemper KE, Xue A, et al. Genotype-by-environment interactions inferred from genetic effects on phenotypic variability in the UK Biobank. Sci Adv. 2019; 5(8): eaaw3538.
25. Zheng H-F, Forgetta V, Hsu Y-H, Estrada K, Rosello-Diez A, Leo PJ, et al. Whole-genome sequencing identifies EN1 as a determinant of bone density and fracture. Nature. 2015; 526(7571): 112-7.
26. Mullin BH, Zhao JH, Brown SJ, Perry JRB, Luan J ’an, Zheng H-F, et al. Genome-wide association study meta-analysis for quantitative ultrasound parameters of bone identifies five novel loci for broadband ultrasound attenuation. Hum Mol Genet. 2017; 26(14): 2791-802.
27. Nabbi A, Almami A, Thakur S, Suzuki K, Boland D, Bismar TA, et al. ING3 protein expression profiling in normal human tissues suggest its role in cellular growth and self-renewal. Vol. 94, European Journal of Cell Biology. 2015. p. 214-22.
28. Chesi A, Wagley Y, Johnson ME, Manduchi E, Su C, Lu S, et al. Genome-scale Capture C promoter interactions implicate effector genes at GWAS loci for bone mineral density. Nat Commun. 2019 19; 10(1): 1260.
29. Maynard RD, Godfrey DA, Medina-Gomez C, Ackert-Bicknell CL. Characterization of expression and alternative splicing of the gene cadherin-like and PC esterase domain containing 1 (Cped1). Gene. 2018 20; 674: 127-33.
30. Määttä JA, Bendre A, Laanti M, Büki KG, Rantakari P, Tervola P, et al. Fam3c modulates osteogenic cell differentiation and affects bone volume and cortical bone mineral density. Bonekey Rep. 2016; 5: 787.
31. Bendre A, Büki KG, Määttä JA. Fam3c modulates osteogenic differentiation by down-regulating Runx2. Differentiation. 2017; 93: 50-7.
32. Martínez-Gil N, Roca-Ayats N, Monistrol-Mula A, García-Giralt N, Díez-Pérez A, Nogués X, et al. Common and rare variants of WNT16, DKK1 and SOST and their relationship with bone mineral density. Sci Rep. 2018 19; 8(1): 10951.
33. Roca-Ayats N, Martínez-Gil N, Cozar M, Gerousi M, Garcia-Giralt N, Ovejero D, et al. Functional characterization of the C7ORF76 genomic region, a prominent GWAS signal for osteoporosis in 7q21.3. Bone. 2019; 123: 39-47.
34. Gonzalez JR, Armengol L, Sole X, Guino E, Mercader JM, Estivill X, et al. SNPassoc: an R package to perform whole genome association studies. Vol. 23, Bioinformatics. 2007. p. 654-5.
35. Mitek T, Nagraba Ł, Deszczyński J, Stolarczyk M, Kuchar E, Stolarczyk A. Genetic Predisposition for Osteoporosis and Fractures in Postmenopausal Women. Adv Exp Med Biol. 2019; 1211: 17-24.
36. García-Ibarbia C, Pérez-Núñez MI, Olmos JM, Valero C, Pérez-Aguilar MD, Hernández JL, et al. Missense polymorphisms of the WNT16 gene are associated with bone mass, hip geometry and fractures. Osteoporos Int. 2013; 24(9): 2449-54.
37. Correa-Rodríguez M, Schmidt Rio-Valle J, Rueda-Medina B. Polymorphisms of the WNT16 gene are associated with the heel ultrasound parameter in young adults. Osteoporos Int. 2016; 27(3): 1057-61.
38. Hendrickx G, Boudin E, Fijałkowski I, Nielsen TL, Andersen M, Brixen K, et al. Variation in the Kozak sequence of WNT16 results in an increased translation and is associated with osteoporosis related parameters. Bone. 2014; 59:57-65.