Genetic variation in WNT16 and its association with bone mineral density, fractures and osteoporosis in children with bone fragility.

Several genome-wide association studies (GWAS), GWAS meta-analyses, and mouse studies have demonstrated that wingless-related integration site 16 (WNT16) gene is associated with bone mineral density (BMD), cortical bone thickness, bone strength and fracture risk. Practically no data exist regarding the significance of WNT16 in childhood-onset osteoporosis and related fractures. We hypothesized that pathogenic variants and genetic variations in WNT16 could explain skeletal fragility in affected children. We screened the WNT16 gene by Sanger sequencing in three pediatric cohorts: 35 with primary osteoporosis, 59 with multiple fractures, and in 95 healthy controls. Altogether, we identified 12 variants in WNT16. Of them one was a rare 5'UTR variant rs1386898215 in genome aggregate and medical trans-omic databases (GnomAD, TOPMED; minor allele frequency (MAF) 0.00 and 0.000008, respectively). One variant rs1554366753, overrepresented in children with osteoporosis (MAF = 0.06 vs healthy controls MAF = 0.01), was significantly associated with lower BMD. This variant was found associated with increased WNT16 gene expression at mRNA level in fibroblast cultures. None of the other identified variants were rare (MAF < 0.001) or deemed pathogenic by predictor programs. WNT16 may play a role in childhood osteoporosis but genetic WNT16 variation is not a common cause of skeletal fragility in childhood.

Introduction

Osteoporosis is a systemic bone disease in which an imbalance between bone resorption and formation results in compromised bone microarchitecture and consequent susceptibility to low-energy fractures. Osteoporosis affects hundreds of millions of people globally and has emerged as a serious health problem, fractures and their complications having devastating effects on an individual and societal level. While osteoporosis typically concerns postmenopausal women, there is growing awareness that it may also affect men and children (Mäkitie, 2013). In fact, fractures are common in children and young adults; studies have found that more than 40% of boys and 25% of girls sustain one or more fractures before the age of 16 (Mäyränpää et al., 2012). We previously conducted a prospective population-based study on fractures in children less than 16 years of age. During a 12-month period 1396 fractures were recorded, of which 63% occurred in boys (Mäyränpää et al., 2010). Most of the fractures in children are low-trauma fractures and occur at the age of 10–14 years (Wasserman and Gordon, 2017).

Several genes are involved in bone mass acquisition, and various monogenic disorders characterized by reduced bone mineral density (BMD) and increased bone fragility have been described (Arundel and Bishop, 2015; Hendrickx et al., 2015; Boudin et al., 2016). Numerous genome-wide association studies (GWAS) have demonstrated that genetic variations in over 30 loci are associated with BMD in children (Zhu et al., 2021).

Although the genetic background and pathogenesis of childhood osteoporosis are still incompletely understood, diligent research has revealed several signal pathways predisposing to skeletal fragility. One of the pathways known to be vital for normal bone metabolism is the Wingless related integration site (WNT) signaling pathway, several of its components have been associated with increased risk of osteoporosis and fractures (Lerner and Ohlsson, 2015; Rivadeneira and Mäkitie, 2016; Baron and Kneissel, 2013).

WNT16 is a key ligand for the pathway and acts mainly on the osteoblastic lineage by regulating both directly and indirectly the maturation and action of osteoclasts (Movérare-Skrtic et al., 2014). Osteoblast expressed WNT16 signals via the canonical WNT pathway to regulate osteoprotegerin (OPG) expression, which in turn functions as a decoy receptor for Receptor activator of nuclear factor-kappa-Β ligand (RANKL) expressed by osteoblasts and osteoclasts and regulates osteoclastogenesis (Gori et al., 2015). The lack of β-catenin in osteoclast precursor cells in mice has been observed to reduce the number and activity of osteoclasts, resulting in osteopetrosis, indicating a direct effect of canonical WNT signaling on osteoclastogenesis (Movérare-Skrtic et al., 2014). The relation between WNT16 and bone homeostasis has been further underlined in several GWAS and GWAS meta-analyses in the over 20 large-scale collaborative studies performed to date. For example, Medina-Gomez et al reported WNT16 variants to associate with whole body BMD in children, Estrada et al identified variants to associate with lumbar spine and femoral BMD, and Zheng et al demonstrated an association between variants and forearm BMD and fractures as well as bone cortical thickness (Medina-Gomez et al., 2018; Estrada et al., 2012; Zheng et al., 2012). These are supported by functional studies in mice, indicating that lack of Wnt16 leads to decreased BMD and an increased risk of non-vertebral fractures (Movérare-Skrtic et al., 2014; Gori et al., 2015). On the other hand, overexpression of Wnt16 does not prevent glucocorticoid treatment induced cortical bone loss in mice (Alam et al., 2018). WNT16 localize near cadherin-like and PC-esterase domain containing 1 (CPED1) and Inhibitors of Growth 3 (ING3) in genome; all three genes have been reported associated with BMD in GWAS and promoter interactions studies (Chesi et al., 2019). So far, no monogenic forms of skeletal fragility resulting from WNT16 mutations have been reported.

Acknowledging the central role of WNT16 in skeletal health and pathogenesis of osteoporosis, we hypothesized that, similar to other components of the WNT signaling pathway such as WNT1 (Laine et al., 2013; Mäkitie et al., 2016) and LRP5 (Saarinen et al., 2010; Gong et al., 2001; Laine et al., 2011), pathogenic variants in WNT16 could play a role in childhood-onset skeletal fragility. The aim of this study was to elucidate the role of WNT16 in the pathogenesis of childhood fractures and childhood-onset osteoporosis by screening for WNT16 variants in three pediatric cohorts: children with primary osteoporosis, children with multiple fractures but normal or decreased BMD, and healthy controls.

Material and methods

Subjects

This study was part of our research program on genetic determinants of metabolic bone disease. The Research Ethics Board at Helsinki University Hospital approved the study, and all patients and/or their legal guardians gave a written consent before study participation. The study included three cohorts: children with primary osteoporosis, children with recurrent fractures, and healthy controls, as outlined below (Table 1). Subjects of the all cohorts were recruited from the same geographical area. Clinical data, including fracture history, were collected by patient interviews and from prior hospital records and imaging studies, as previously described (Saarinen et al., 2010; Mäyränpää et al., 2011; Pekkinen et al., 2012). Fractures caused by impacts comparable to a fall of ≥3 m were considered high-energy fractures.

Table 1

Clinical characteristics of the study cohorts.

Characteristics	Osteoporosis patients	Fracture patients	Controls (healthy school children)
Number of subjects	35	59	95
Age, median (range; years)	12.0 (6.7–16.6)	10.2 (4.4–16.8)	11.7 (7–19)
Males	20 (57%)	40 (69%)	45 (48%)
Inclusion criteria	(i) Clinically diagnosed primary osteoporosis, (ii) No known genetic cause of OI or primary osteoporosis, (iii) Secondary causes of osteoporosis excluded, (iv) Low BMD (Z-score below −2.0) and/or at least 3 peripheral fractures and/or at least 1 vertebral compression fracture caused by low- or moderate-energy trauma	i) Age 4–16 years, ii) A clinically significant fracture history with ≥2 low-energy long-bone fractures before age 10 years, ≥3 low-energy long-bone fractures before age 16 years, or ≥1 low-energy vertebral fracture (loss of ≥20% vertebral height), iii) No diagnosis or suspicion of OI and no underlying disease likely to explain their bone fragility	Healthy school children without any underlying severe chronic illness.
Subjects with moderate decline in BMD (BMD Z-score between −1.0 and −2.0)
Lumbar spine	9 (26%)	12 (20%)	11 (12%)
Femoral	5 (14%)	5 (9%)	8 (8%)
Whole body	7 (20%)	7 (12%)	7 (7%)
Subjects with severe decline in BMD (BMD Z-score < −2.0)
Lumbar spine	13 (37%)	4 (7%)	1 (1%)
Femoral	20 (59%)	2 (3%)	0 (0%)
Whole body	5 (23%)	1 (2%)	0 (0%)
Subjects with low- or moderate-energy trauma peripheral fractures	24 (69%)	54 (92%)	18 (19%)
Subjects with compression fractures	21 (60%)	11 (19%)	0 (0%)
Publication of original study	(Mäyränpää et al., 2011)	(Mäyränpää et al., 2012)	(Pekkinen et al., 2012)

BMD = bone mineral density.

Cohort 1: Osteoporosis cohort

The cohort consisted of 35 patients (20 boys and 14 girls; age range 4–17 years, median 12 years) recruited at the Metabolic Bone Clinic, Children's Hospital, Helsinki University Hospital, Finland over a 10-year period. The inclusion criteria were (1) clinically diagnosed primary osteoporosis, (2) exclusion of type I collagen-related osteogenesis imperfecta (OI), either clinically or by genetic testing, and (3) exclusion of secondary osteoporosis with biochemical and individually determined clinical evaluations. Further, all had been screened and found negative for pathogenic variants in WNT1 and LRP5, and several had undergone more extensive genetic testing before inclusion. Diagnosis of osteoporosis was based on BMD measurements (Z-score below −2.0) and/or a clinically significant fracture history with at least three peripheral fractures and/or at least one vertebral compression fracture caused by low- or moderate-energy trauma (Mäkitie, 2013; NIH Consensus Development Panel on Osteoporosis Prevention, 2001). Bone biopsies were obtained from several patients to confirm the diagnosis of osteoporosis. Biopsy was done according to standard methods: in brief, following a double labelling course with oral tetracycline (2 days plus 2 days with a 10-day interval), biopsies were obtained, in local anesthesia, from the anterior iliac spine with a manual drill with a trephine of 7.5 mm inner diameter (Rochester Bone Biopsy; Medical Innovations International, Rochester, MN). All had clinically diagnosed childhood-onset primary osteoporosis. None of the patients had a family history or a known diagnosis of OI, or a genetically identified molecular cause of their disease.

Cohort 2: Fracture cohort

The cohort consisted of 59 children with multiple fractures. All were recruited at the Children's Hospital, Helsinki, Finland, during a prospective epidemiological study (Mäyränpää et al., 2012; Mäyränpää et al., 2010). Over 12 months, all children (N = 1412) aged 4–16 years who had been treated for an acute and radiographically confirmed fracture before the age of 16 were assessed for their fracture history. Trauma mechanisms and previous medical histories were available for 1361 (96%) of the children. The inclusion criteria for Fracture cohort were (1) age 4–16 years, (2) ≥2 low-energy long-bone fractures before age 10 years, (3) ≥3 low-energy long-bone fractures before the age of 16 years, or (4) ≥1 low-energy vertebral fracture (loss of ≥20% vertebral height). Children with a diagnosis or suspicion of OI were excluded, as were children with an underlying disease likely to explain their bone fragility.

From the 1412 children, altogether 71 patients fulfilled these criteria; DNA from peripheral blood was available for 64 (Kämpe et al., 2017). Of them, five fulfilled the criteria for primary osteoporosis and these patients were subsequently categorized into the Osteoporosis cohort. Thus, the Fracture cohort included 59 children with 19 girls and 40 boys, ranging in age from 4 to 16 years (median 10.2 years) (Table 1). Data on blood biochemistry, spinal radiographs, and DXA measurements were collected for all participants.

Cohort 3: Control cohort

The healthy controls included 95 healthy children (45 boys and 50 girls) aged 7–19 years (median 11.7 years), without any evident skeletal disorder or other underlying severe chronic illnesses (Table 1). They were initially recruited to a school-based study investigating vitamin D status and its association with bone health in Finnish school children (Pekkinen et al., 2012).

Dual-energy X-ray absorptiometry measurements

BMD, bone mineral content (BMC), and bone area (BA) were measured with DXA (Hologic Discovery A, pediatric software, version 12.4) from the lumbar spine (LS; L1–L4), total hip (FEM) and whole body (WB). All measured values were transformed into Z-scores using equipment-specific age- and sex-adjusted reference data for Caucasian US children (Kelly et al., 2009); all subjects were of normal height for age.

DNA extraction and sequencing

DNA was extracted from whole blood samples using standard procedures. The WNT16 gene (NC_00000007.14), including all five exons (exon 1A, 1B, 2, 3, 4), a minimum of 30 bases of flanking introns, and 5′ and 3′ untranslated (UTR) regions, was Sanger sequenced in all subjects. Primers's sequences are shown in Table 2S. Chromatograms were analyzed with Sequencher 5.0 (Gene Codes Corporation, Ann Arbor, MI, USA). Information regarding detected single nucleotide variants (SNPs) were collected from Ensembl (https://www.ensembl.org/), NCBI databases (https://www.ncbi.nlm.nih.gov/) (Ensembl – Home Page, n.d.; NCBI, n.d.) (Table 2), and predictions performed by Varsome (Kopanos et al., 2019) with the following prediction programs: DANN, DEOGEN2, EIGEN, FATHMM-MKL, M-CAP, MVP, MutationAssessor, MutationTaster, PrimateAI, REVEL, SIFT and GERP.

Table 2

Genetic variation in WNT16 in the study populations.

SNP	Change (NM_016087.2)	Variant type	Position in gene	MAF in cohorts	MAF in gnomAD ALL	MAF in gnomAD Finns	ACMG classification and prediction	Previous reports of association with skeletal health
rs1554366753	c.-6 -5insTACC	Untranslated	5′UTR/Exon 1A	0.06a 0.0b 0.01c	0.008406	0.01373	Likely benign
rs10668066	c.-3_1insCACC, (p.Met1ThrfsTer141)	Untranslated	5′UTR/Exon1A	0.3a 0.21b 0.025c	0.2348	0.2178	Benign	(Kemp et al., 2017; Kim, 2018)
rs1386898215	c.-246T > A	Untranslated	5′UTR/Exon 1B	0.01a 0.0b 0.0c	0.0	0.0	Uncertain significance
rs201022838	c.-15C > A	Untranslated	5′UTR/Exon 1B	0.0a0.009b 0.0c	0.001134	0.0002972	Likely benign	(Martínez-Gil et al., 2018)
rs35391640	c.66C > G, (p.Leu22Leu)	Synonymous	Exon 1B	0.07a 0.07b 0.07c	0.02980	0.07725	Benign	(Martínez-Gil et al., 2018)
rs2908004	c.244G > A, (p.Gly82Arg)	Missense	Exon 2	0.44a 0.37b 0.42c	0.4567	0.4188	Benign	(Kim, 2018; Kemp et al., 2014)
rs17143291	c.300C > A, (p.Thr100Thr)	Synonymous	Exon 2	0.04a 0.0b 0.021c	0.009273	0.01395	Likely benign	(Martínez-Gil et al., 2018)
rs17143296	c.630G > A, (p.Arg210Arg)	Synonymous	Exon 3	0.0a 0.03b 0.011c	0.009818	0.004900	Likely benign	(Martínez-Gil et al., 2018)
rs74389152	c.631C > A, (p.Gln211Lys)	Missense	Exon 3	0.0a 0.008b 0.02c	0.001923	0.006333	Likely benign
rs2707466	c.788C > T, (p.Thr263Ile)	Missense	Exon 4	0.44a 0.37b 0.04c	0.4550	0.4134	Benign	(Martínez-Gil et al., 2018)
rs17143305	c.*113C > T	Untranslated	3′UTR/Exon4	0.08a 0.05b 0.01c	0.1374	0.1550	Benign	(Martínez-Gil et al., 2018)
rs3832519	c.*128_*129insCTCT	Untranslated	3′UTR/Exon4	0.13a 0.20b 0.15c	0.2588	0.1663	Benign	(Martínez-Gil et al., 2018)

MAF = minor allele frequency.

NC_00000007.14 = WNT16 gene reference sequence.

ACMG = The American College of Medical Genetics and Genomics.

Deviating MAF values in cohorts are bolded.

Primary osteoporosis patients.

Fracture patients.

Healthy controls.

RNA extraction and gene expression analyses

Fibroblast cells for RNA extraction were available for Patient 1 and three age and gender-matched controls. Cells were cultured from skin biopsies according to our previous protocol (Vakkilainen et al., 2019). RNA was isolated using the RNAeasy Mini kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. RNA concentrations were measured using a Qubit 2.0 Fluorometer (Thermo Fisher, Carlsbad, CA). cDNA was transcribed from 1 μg of total RNA using the QuantiTect Reverse Transcription Kit cDNA synthesis kit (Qiagen) according to the manufacturer's protocol. qRT-PCR assays were performed in quintuplicates using the CXF96 Real-Time system (Bio-Rad Laboratories, Hercules, CA) and using the TaqMan WNT16, Actin Beta (ACTB) and TATA-Box Binding Protein (TBP) gene expression assays Hs00365138_m1, Hs01553299_g1, Hs99999903_m1 and Hs00427620_m1, respectively (Applied Biosystems, Thermo Fisher Scientific, Carlsbad, CA). ACTB and TBP were used as reference genes for data normalization. Threshold cycle (Ct) values were determined using CFX Manager Software (Bio-Rad Laboratories). Relative expression was calculated using the comparative Ct or 2ΔΔCt method (Livak and Schmittgen, 2001; Pfaffl, 2001).

Statistical methods

All variables were checked and tested for normality. Differences between groups were tested with Independent Samples t-test or Mann-Whitney U test, as appropriate. Categorical variables were tested with chi-square test. Analysis of variance (ANOVA) or a Kruskal-Wallis test was applied to compare more than two groups, posthoc comparisons were made with Bonferroni correction.

Analysis of covariance (ANCOVA) was used in examining associations between SNPs and BMD or fractures in all cohorts. Association between variables was tested with Pearson or Spearman correlation. All tests were two-tailed, and p-values <0.05 were considered significant. Calculations were performed using SPSS Statistics 25 (IBM Corporation). ANOVA, a Mann-Whitney U test or a Kruskal-Wallis test was used to test the association between genetic variation and LS, FEM, and WB BMD Z-scores.

Results

Overall genetic variations in the WNT16 gene

Altogether, in all three cohorts (Osteoporosis cohort N = 35, Fracture cohort N = 59, Control cohort N = 95), we identified 12 variants in the WNT16 gene (Table 2, Table 1S). Of them, six were present in all cohorts, one only in the Osteoporosis cohort, one only in the Fracture cohort, and none only in the Control cohort. The five that were detected in homozygous state are also present in heterozygous state. Six of the variants were located in the coding region (two in exon 2), and six in untranslated regions (Fig. 1).

Fig. 1

Locations of the 12 genetic changes found in WNT16 in the three pediatric cohorts. Locations of the 12 genetic changes found in WNT16. Untranslated regions are shown in light blue and translated regions in dark blue. Synonymous = not cause amino acid change.

None of the 12 variants were novel. However, one of them was a rare untranslated variant rs1386898215 (gnomAD, MAF = 0.00 and TOPMED MAF = 0.000008), located in the 5′UTR region and found heterozygous in one of the Osteoporosis patients (1/35, 3%). The reference SNP (rs) numbers for the other identified variants were rs1554366753, rs10668066, rs201022838, rs35391640, rs2908004, rs17143291, rs17143296, rs74389152, rs2707466, rs3832519 and rs17143305 (Table 2). Of these, three were synonymous SNPs, three nonsynonymous SNPs, and six in the untranslated regions (Fig. 1). None of the variants were considered disease-causing according to prediction programs. Details and allele frequencies for all identified variants are presented in Tables 2, 1S.

WNT16 genetic variation in the Osteoporosis cohort and the Fracture cohort

For some of the variants, the allele frequencies varied between the groups. Rs1554366753 minor allele was relatively more common in the Osteoporosis cohort (6%, vs 0% in the Fracture cohort and 1% in the Control cohort) and associated with primary osteoporosis (khi2, p = 0.048). The variation is a four-nucleotide insertion in the immediate proximity of exon 1A. The first nucleotide varied; while the change insCACC, with identifier number rs10668066, was common in all cohorts, insTACC (rs1554366753) was uncommon and present heterozygous only in four primary osteoporosis patients (Osteoporosis cohort) and two healthy children (Control cohort). Based on the higher frequencies of some of the variants in the Osteoporosis and Fracture cohorts, compared with the Control cohort, MAF in the global population or the Finnish population, three changes (rs1554366753/rs10668066, rs1386898215, and rs201022838) were considered significant. All these variants were in the 5′UTR region of WNT16. The clinical features of the patients with these WNT16 variants are shown in Table 3. Patients 1, 2, 3 and 4 are children with primary osteoporosis (Osteoporosis cohort), and patient 5 is a fracture-prone child (Fracture cohort). All patients have had peripheral fractures, and all the Osteoporosis patients had BMD Z-score values below −2.0 before bisphosphonate medication.

Table 3

Clinical features of patients with significant genetic changes in WNT16.

Patient	Cohort	SNP	Variant location	Age at last follow up	Sex	BMD prior to BP (Z-score)			Number of peripheral fractures	Compression fractures
						LS	FEM	WB
1	Osteoporosis	rs1554366753	5′UTR	19	M	−2.9	−2.2	−2.8	6	Yes
2	Osteoporosis	rs1554366753	5′UTR	13	F	−1.8	−1.8	−0.5	5	No
3	Osteoporosis	rs1554366753	5′UTR	13.6	M	−1.7	−1.8	−1.6	1	Yes
4	Osteoporosis	rs1386898215	5′UTR	20	M	−4.1	−3.3	–	4	Yes
5	Fracture	rs201022838	5′UTR	6.5	M	−0.1	0.7	1.1	7	No

BMD values have been measured before patients have started bisphosphonate (BP) medication.

BMD = bone mineral density.

Genetic variation in WNT16 and its association with bone mineral density and fractures

The SNP rs1554366753 was associated with BMD. The minor allele insTACC was linked to lower BMD in FEM in the Osteoporosis cohort compared with allele frequencies and FEM Z-score values in the Control cohort (Kruskal-Wallis, p = 0.038). In contrast, there were no statistically significant differences between LS and WB Z-scores or peripheral fractures and genetic variation in rs1554366753. The rs10668066 genotypes and MAF frequencies differed between the Osteoporosis and Control cohorts (Table 2, Kruskal-Wallis, p = 0.036). However, no statistically significant associations between rs10668066 and BMD or fractures were observed (ANCOVA, LS p = 0.666, FEM p = 0.600 and WB p = 0.343, peripheral fractures p = 0.723) using cohorts as covariate. Furthermore, the other ten studied WNT16 SNPs did not associate with BMD (FEM, LS and WB Z-scores), osteoporosis or peripheral fractures (Table 4).

Table 4

Association of WNT16 genetic variation between bone mineral density and peripheral fractures in all study cohorts (N = 189).

SNP	LS BMD Z score (mean ± SD)				Femoral BMD Z-core (mean ± SD)				Whole body BMD (mean ± SD)				Peripheral fractures (mean ± SD)
	AA	Aa	aa	P	AA	Aa	aa	P	AA	Aa	aa	P	AA	Aa	aa	P
rs1554366753	−0.43 ± 1.25	−1.55 ± 1.29		0.120a	−0.62 ± 1.06	−1.38 ± 1.13		0.240a	−0.08 ± 0.90	−1.13 ± 1.38		0.160a	1.36 ± 1.73	3.00 ± 2.94		0.211a
rs10668066	−0.44 ± 1.23	−0.40 ± 1.33	−1.46 ± 1.13	0.666b	−0.11 ± 1.03	−0.01 ± 1.13	−0.96 ± 1.35	0.600b	−0.90 ± 0.95	−0.28 ± 0.83	−1.02 ± 1.22	0.343b	1.42 ± 1.57	1.26 ± 2.00	2.60 ± 2.70	0.723b
rs1386898215	−0.43 ± 1.23	−1.40 ± 0.00		0.120a	−0.07 ± 1.05	−3.3 ± 0.00		0.240a	−0.10 ± 0.93			0.155a	1.38 ± 1.76	4.00 ± 0.00
rs201022838	−0.45 ± 1.26	−0.10 ± 0.00		0.831a	−0.95 ± 1.08	−0.70 ± 0.00		0.451a	−0.11 ± 0.93	1.10 ± 0.00		0.197a	1.36 ± 1.72	7.00 ± 0.00
rs35391640	−0.42 ± 1.25	−0.70 ± 1.37	0.05 ± 0.35	0.864b	−0.12 ± 1.09	0.06 ± 1.04	0.55 ± 0.21	0.385b	−0.10 ± 0.91	−0.13 ± 1.13	0.05 ± 0.07	0.933b	1.31 ± 1.60	2.14 ± 2.62	0.00 ± 0.00	0.167b
rs2908004	−0.50 ± 1.14	−0.40 ± 1.24	−0.64 ± 1.58	0.850b	−0.17 ± 1.01	−0.02 ± 1.01	−0.18 ± 1.09	0.585b	−0.13 ± 0.90	−0.06 ± 0.90	−0.18 ± 1.09	0.848b	1.37 ± 1.60	1.38 ± 1.90	1.5 ± 1.69	0.928b
rs17143291	−0.44 ± 1.28	−0.82 ± 1.01		0.351a	−0.80 ± 1.08	−0.52 ± 1.20		0.480a	−0.09 ± 0.93	−0.34 ± 088		0.586a	1.40 ± 1.76	1.40 ± 2.07		0.964a
rs17143296	−0.46 ± 1.27	0.20 ± 0.72		0.244a	−0.11 ± 1.08	0.56 ± 1.12		0.284a	−0.11 ± 0.93	0.35 ± 0.58		0.312a	1.39 ± 178	1.5 ± 1.00		0.603a
rs74389152	−0.45 ± 1.27	−0.27 ± 0.91		0.796a	−0.09 ± 1.08	0.03 ± 1.01		0.952a	−0.10 ± 0.93	−0.00 ± 0.95		0.991a	1.40 ± 1.76	1.33 ± 2.31		0.805a
rs2707466	−0.46 ± 1.12	−0.41 ± 1.24	−0.57 ± 1.62	0.956b	−0.18 ± 0.99	−0.01 ± 1.02	−0.15 ± 1.41	0.502b	−0.13 ± 0.89	−0.06 ± 0.91	−0.14 ± 1.08	0.852b	1.32 ± 1.53	1.36 ± 1.91	1.69 ± 1.78	0.953b
rs17143305	−0.47 ± 1.30		0.10 ± 0.46	0.230a	−0.11 ± 1.08		0.60 ± 0.80	0.140a	−0.12 ± 0.93		0.52 ± 0.44	0.072a	1.35 ± 1.74		3.0 ± 2.00	0.067a
rs3832519	−0.48 ± 1.17	−0.42 ± 1.52	0.10 ± 0.46	0.237b	−0.13 ± 1.06	−0.07 ± 1.14	−0.60 ± 0.80	0.093b	−0.19 ± 0.95	−0.07 ± 0.87	0.52 ± 0.44	0.022b	1.32 ± 1.81	1.42 ± 1.58	3.0 ± 2.00	0.165b

p = p-values.

BMD = bone mineral density.

LS = lumbar spine, FEM = femoral, WB = whole body.

AA = major allele homozygotes, Aa = heterozygous, aa = minor allele homozygous.

Kruskal Wallis test.

ANCOVA, cohort as a covariant.

WNT16 gene expression

In order to analyze whether the rs1554366753 minor allele variant affects gene expression, we compared WNT16 gene expression in Patient 1, harboring the heterozygous variant insTACC, and three healthy age- and gender-matched controls without the variant. WNT16 mRNA expression levels in cultured fibroblasts were four to five times higher in the patient than in the controls in both WNT16 Taqman assays (Fig. 2).

Fig. 2

Normalized WNT16 mRNA expression measured with two Taqman assays (WNT16_g and WNT16_m) in Osteoporosis patient 1, heterozygous for the WNT16 variant insTACC (rs1554366753), vs. three age- and sex-matched healthy controls. *** P < 0.001.

Discussion

The key genetic factors relating to childhood osteoporosis remain poorly understood. Several GWAS and GWAS meta-analyses, mainly in older adults, have shown that WNT16 is associated with BMD, cortical bone thickness, and fracture risk (Zheng et al., 2012; Kämpe et al., 2017; Ensembl – Home Page, n.d.; Liu et al., 2020). Extensive studies on Wnt16 null mice have further demonstrated that WNT16 is a vital regulator of cortical bone mass and strength (Lerner and Ohlsson, 2015; Movérare-Skrtic et al., 2014; Zheng et al., 2012). Data on the significance of WNT16 in childhood-onset osteoporosis and risk of fractures is, however, still largely lacking.

Here, we describe, for the first time, the prevalence and nature of WNT16 variants in pediatric patients. We screened the WNT16 gene by Sanger sequencing in 35 pediatric patients with primary osteoporosis and 59 pediatric patients with recurrent fractures and hypothesized that pathogenic variations in WNT16 would explain the clinical characteristics in of some of these patients. In total, we found 12 variants, of which none were regarded pathogenic and most were known and predicted likely benign. Overall, the allele frequencies were similar between patient cohorts and healthy controls and comparable to database frequencies, suggesting that variants in WNT16 are not a common cause of pediatric osteoporosis or fragility fractures.

Despite this, three interesting variants were found. The first of them was a heterozygous change c.-6 -5insTACC (rs1554366753) in the 5′UTR region of WNT16. InsTACC had a significantly higher allele frequency in Osteoporosis patients (p = 0.038) than in healthy children, the overall population, or the Finnish population (Ensembl – Home Page, n.d.). Furthermore, mRNA expression of this variant was elevated in heterozygous Patient 1 compared to Healthy controls. In addition, insCACC allele frequencies also differed between the Osteoporosis patients and Healthy controls (p = 0.036). Previously, this variant has been reported to associate with heel bone mineral density in GWAS studies (Kemp et al., 2017; Kim, 2018). A transcriptome-wide association study recently detected 88 genes significantly associated with total body BMD or fracture through expression or ribonucleic acid splicing; WNT16 was one of the reported genes (Liu et al., 2020).

Secondly, we identified a rare change c.-246T > A (rs1386898215) in the WNT16 5′UTR region in one Osteoporosis patient (Patient 4). The heterozygous change was absent in the other two cohorts (Table 1S) and has not been reported in any other database, report or GWAS study. Since the change was present in only one patient and was in the 5′UTR region, it is likely uncommon in patients with primary osteoporosis.

Lastly, the variant rs201022838 was present in heterozygous only one patient in the Fracture cohort (Patient 5). This c.-15C > A change is in the 5′UTR region, and the allele frequency was higher in this cohort than in the other two cohorts, the overall population, or the Finnish population. Although frequency difference was not statistically significant, some of this may be accounted for by the small cohort sizes.

Some of the identified variants have previously been reported in other studies in older subjects. Ibarbia et al found that variants rs2707466 in exon 4 and rs2908004 in exon 2 strongly associate with BMD and fracture risk in individuals under the age of 80 (García-Ibarbia et al., 2013). Similarly, Correa-Rodríquez et al found that variant rs2908004 is associated with broadband ultrasound attenuation (BUA) measurement and concluded that, in young adults, WNT16 affects bone mass (Correa-Rodríguez et al., 2016). Another study showed that variant rs2707466 affects BMD at several sites. However, they found rs2908004 to be only mildly associated with BMD (Hendrickx et al., 2014). Zheng et al found that rs2707466 had genome-wide significant association with forearm BMD and cortical thickness in females and males under 20 years of age (Zheng et al., 2012). Lastly, Kemp et al have shown an association between rs2908004 and upper limb BMD in a pediatric study population (Kemp et al., 2014). In our study, rs2707466 and rs2908004 were present at similar frequencies in all three cohorts with MAF 0.44, 0.37 and 0.40, respectively. Therefore, they are not associated with BMD in our pediatric patients. Hendrickx et al also found that rs55710688 affects translation efficiency and the amount of Wnt16 and, therefore, BMD. This study included over 1000 individuals, bringing more statistical power to the analyses. These results are in line with our findings. We found rs55710688, whose allelic variations are currently part of the rs10668066 locus, to be present at different frequencies in the Osteoporosis patients and Healthy controls. Furthermore, rs10668066 has been associated with heel bone mineral density in adults (Kim, 2018; Kemp et al., 2014). In a recent GWAS study, Surakka et al reported that the intergenic WNT16 variant rs2707518 was one of ten genome-wide significant BMD-associated loci in the discovery HUNT dataset (N = 19,705) (Surakka et al., 2020).

Wnt16 null mice (Wnt16−/−) have increased cortical porosity, reduced cortical bone trajectory and bone strength; therefore, there is greater fracture susceptibility at cortical sites (Zheng et al., 2012). WNT16 is mainly expressed in osteoblasts and affects osteoclastogenesis both directly and indirectly (Movérare-Skrtic et al., 2014). WNT16 exerts its effects on bone through osteoblastic lineage cells and is an agonist of both the canonical and non-canonical pathways (Movérare-Skrtic et al., 2014). Wnt16 mRNA expression is most significant at endosteal surfaces of cortical bone and proper cortical bone homeostasis is dependent on the expression of Wnt16 by early osteoblasts during both development and skeletal growth. Contrary to this, Wnt16 does not affect trabecular bone homeostasis (Movérare-Skrtic et al., 2014; Gori et al., 2015).

As seen in Wnt16−/− mice, the number of osteoblasts remains unaltered, but osteoclastogenesis is increased at endocortical surfaces (Movérare-Skrtic et al., 2014). Wnt16 has an anabolic effect on bone. In our study, the mRNA expression of WNT16 was elevated in Patient 1, who has a relatively rare rs1554366753 (MAF < 0.015 in Finns and 0.0008 in gnomAD ALL) insTACC variant in the 5′UTR region, compared to sex- and age-match controls. This result is opposite to that seen in the knockout mice. However, we hypothesize that constantly high WNT16 expression levels can affect bone homeostasis, especially the endothelial surface of the cortical bone, where WNT16 mRNA expression is high. Previous studies have noticed that a certain level of Wnt16 is required for bone homeostasis (Zheng et al., 2012; García-Ibarbia et al., 2013), while Wnt16 overexpression might inhibit osteoblast differentiation (Jiang et al., 2014).

Patient 1 had low (<−2.0) BMD in LS, FEM, and WB and many peripheral fractures. It remains to be explored in future studies and bigger sample sizes whether the relatively rareWNT16 variant rs1554366753 and WNT16 overexpression contribute to skeletal fragility.

We recognize several limitations in the study. The patient cohorts screened were rather small; a larger number of patients could have resulted in more genetic findings and enabled us to perform more detailed statistical analyses. Also, this study only screened WNT16 exons and UTR regions; regulatory sequences and introns remain to be screened. However, this study was not designed to only detect associations but to also screen for rare disease-causing variants. Furthermore, not all patients screened had a severe illness, and pathogenic WNT16 variants could be perhaps discovered in patients with more severe osteoporotic findings. On the other hand, to the best of our knowledge, no known WNT16 pathogenic variants causing osteoporosis have been found. Despite these limitations, this is the largest group of pediatric patients screened for genetic changes in WNT16. Therefore, the results are significant and novel. These findings need to be confirmed in other pediatric Osteoporosis and Fracture cohorts in the future. It may be that WNT16 variants play a more significant role in the development of osteoporosis in older patients due to the cumulative effects of aging. Similar studies on larger cohorts are encouraged to further assess the prevalence of these genetic changes and their effect on bone health on the population level.

Conclusion

Our findings indicate that pathogenic variants in WNT16 are not a common cause of increased fractures or osteoporosis in children. Taking into consideration recent contrasting observations demonstrating that the WNT16 genetic variation impacts cortical bone thickness, BMD, bone strength, and fracture risk in mice and humans, it may be that rare WNT16 variants do not play an essential role in the pathogenesis of childhood osteoporosis but serve a more central role later in adulthood. Lastly, as the WNT signaling pathway has become an interesting target for new osteoporosis medications, comprehensive knowledge on the functions of WNT16 in bone health and disease in all phases of life is valuable for improved diagnostic and therapeutic means.

The following are the supplementary data related to this article.

Table 1S

WNT16 genotype frequencies in the Osteoporosis and Fracture cohorts and the Healthy controls.

Table 2S

Wnt16 primers.

CRediT authorship contribution statement

Riikka Mäkitie: Data collection, Drafting of the manuscript and approving the final version of the manuscript. Sara Mäkitie: Data collection and Analysis, Drafting of the manuscript and Approving the final version of the manuscript. Mervi Mäyränpää: Data collection, Revising of the manuscript and approving the final version of the manuscript. Minna Pekkinen: Study design, Data collection and analysis, Drafting of the manuscript and approving the final version of the manuscript, Supervision.

Declaration of competing interest

The authors declare that they have no competing interests.