Comprehensive Genetic Analysis of 128 Candidate Genes in a Cohort With Idiopathic, Severe, or Familial Osteoporosis.

CONTEXT: The genetic bases of osteoporosis (OP), a disorder with high heritability, are poorly understood at an individual level. Cases of idiopathic or familial OP have long puzzled clinicians as to whether an actionable genetic cause could be identified.
OBJECTIVE: We performed a genetic analysis of 28 cases of idiopathic, severe, or familial osteoporosis using targeted massively parallel sequencing.
DESIGN: Targeted sequencing of 128 candidate genes was performed using Illumina NextSeq. Variants of interest were confirmed by Sanger sequencing or SNP array. PATIENTS AND
SETTING: Thirty-seven patients in an academic tertiary hospital participated (54% male; median age, 44 years; 86% with fractures), corresponding to 28 sporadic or familial cases. MAIN OUTCOME MEASURE: The identification of rare stop-gain, indel, splice site, copy-number, or nonsynonymous variants altering protein function.
RESULTS: Altogether, we identified 28 variants of interest, but only 3 were classified as pathogenic or likely pathogenic variants: COL1A2 p.(Arg708Gln), WNT1 p.(Gly169Asp), and IDUA p.(His82Gln). An association of variants in different genes was found in 21% of cases, including a young woman with severe OP bearing WNT1, PLS3, and NOTCH2 variants. Among genes of uncertain significance analyzed, a potential additional line of evidence has arisen for GWAS candidates GPR68 and NBR1, warranting further studies.
CONCLUSIONS: While we hope that continuing efforts to identify genetic predisposition to OP will lead to improved and personalized care in the future, the likelihood of identifying actionable pathogenic variants in intriguing cases of idiopathic or familial osteoporosis is seemingly low.

Osteoporosis (OP) is a common disease characterized by low bone mass and microarchitectural deterioration of bone tissue, resulting in increased fracture risk and hence great morbidity, mortality, and financial burden to health care systems. OP is generally a multifactorial disorder associated with several clinical risk factors such as advanced age, postmenopausal status, low body mass index, and use of glucocorticoids, amongst others [1]. Family history of OP and/or fragility fracture is also a well-recognized clinical risk factor, highlighting the genetic component of bone fragility. In fact, the heritability of bone mineral density (BMD) in the general population has been estimated to be between 50% and 85% [2, 3]. As in other multifactorial traits, the genetic architecture of OP may involve a combination of common allelic variants with low individual phenotypic impact or single, rare, high-impacting variants.

The contribution of common low-impacting genetic variants on BMD has been scrutinized in the last decade by genome-wide association studies (GWAS) [4-29]. In total, these have associated more than 220 loci and 240 genes with low BMD and fractures [29, 30]. Although the biological mechanisms by which these loci contribute to bone fragility remain to be determined in many cases, GWAS have accurately identified well-established factors involved in the regulation of bone mass such as members of the WNT and RANK signaling pathways.

The potential of single, rare, genetic variants to determine striking bone fragility phenotypes is represented by Mendelian disorders such as osteogenesis imperfecta (OI) and Hajdu-Cheney syndrome, caused by high-impacting variants determining severe bone fragility, generally from birth, and often accompanied by extraskeletal features. Monogenic defects in at least 17 genes have been associated with OI, most frequently in COL1A1 or COL1A2, whereas heterozygous variants in NOTCH2 have been shown to cause Hajdu-Cheney syndrome.

In this context, familial forms of OP or severe cases of bone fragility in the absence of recognizable risk factors have long puzzled clinicians. A genetic etiology has been occasionally sought in these cases, leading to the identification of defects in LRP5, WNT1, and PLS3 [31-34]. Interestingly, WNT1 and LRP5 are also involved with more severe bone fragility in OI and OP-pseudoglioma syndrome, respectively, demonstrating that a phenotypic spectrum of allelic variation in these genes exists.

Advances in gene sequencing technology, namely the availability of massively parallel sequencing (MPS), now enable comprehensive analysis of the molecular basis of OP. Several candidate genes can be analyzed simultaneously, allowing not only the inclusion of a broader set of candidates—for example, those arising from GWAS—but also the recognition of oligogenic interactions as disease mechanisms. The aim of this study was to identify genetic variants associated with idiopathic, severe, or familial forms of OP by targeted sequencing of 128 candidate genes, hoping to improve our understanding of the complex genetic architecture of OP, and, on an individual level, to refine the clinical management of patients and their relatives.

Methods

Participants

This study included individuals with idiopathic, severe, and/or familial OP/low BMD. In accordance with the 2015 International Society for Clinical Densitometry (ISCD) Official Positions, OP was defined based on BMD assessment by dual-energy x-ray absorptiometry (DXA) in postmenopausal women and men age 50 years and older by a T score of –2.5 or less on the lumbar spine (LS), femoral neck (FN), or total hip (TH). In premenopausal women and men younger than 50 years, low BMD was defined by a Z score of –2.0 or less. Men or premenopausal women with idiopathic OP or low BMD were included, as well as postmenopausal women with OP if their T score was –4.0 or less or if a first-degree relative was included for other criteria.

We have studied a Brazilian cohort, the majority of which was recruited from a tertiary OP clinic at the Endocrine Division of Hospital das Clinicas HCFMUSP, Faculdade de Medicina, Universidade de São Paulo, in São Paulo, Brazil. Additional cases were referred by endocrinologists and rheumatologists practicing in São Paulo, Brazil.

All individuals underwent a systematic evaluation consisting of a detailed medical history, physical examination, and laboratory assessment with complete blood count, determination of serum levels of calcium, phosphate, parathyroid hormone, 25-hydroxyvitamin D, alkaline phosphatase, albumin, creatinine, liver enzymes, thyrotropin, free thyroxine, C-reactive protein, total immunoglobulin A, antiendomyseal antibody, and testosterone (for men), and 24-hour urinary calcium. Individuals with secondary osteoporosis as identified by this evaluation were excluded from this study.

A detailed history of nonvertebral fractures was obtained from all participants. Vertebral fractures were ascertained through lateral thoracic and lumbar spine radiographs in all but 4 individuals, using the semiquantitative technique. All patients underwent DXA (Hologic or GE Lunar) according to ISCD standards.

All participants provided informed consent; the study was carried out according to the principles of the Declaration of Helsinki and was approved by the local ethics committee (Comissão de Ética para Análise de Projetos de Pesquisa, Hospital das Clinicas HCFMUSP, approval number 294142).

DNA Extraction

DNA was extracted from peripheral blood or saliva samples. For peripheral blood leukocyte DNA extraction, an in-house salting-out method was used. Salivary DNA was obtained using the Oragene-DNA OG-500 kit (DNA Genotek Inc), following the manufacturer’s instructions. All DNA samples were submitted to quality control before further genetic analyses.

Gene Panel Selection

We conducted a systematic search for genes previously associated to bone strength by GWAS, Mendelian disorders, or nonsyndromic idiopathic/familial OP. Our approach to identifying candidate genes for bone fragility has been previously described [30]. Briefly, the descriptors “osteoporosis,” “fractures,” “bone fragility,” “BMD,” “genes,” “genetics,” and “GWAS” were queried in the PubMed and Online Mendelian Inheritance in Man search engines to identify published studies. From GWAS, we selected for inclusion in the panel candidate genes located 500 kb upstream or downstream from significant (P value < 5 × 10–8) risk loci that were either nominally identified in the original studies or presented additional evidence of involvement in bone metabolism arising from human diseases, mouse models, or in vitro studies [4-28]. Candidate genes associated with Mendelian disorders incurring in low BMD were included, as well as those associated with high BMD, based on the premise of potential phenotypic spectra (for example, distinct variants in LRP5 have been associated with either high or low BMD). Finally, genes that had been associated with idiopathic, pregnancy-associated, or familial OP in candidate gene or whole-exome sequencing studies were also included. Following this systematic search, in March 2015 a panel was designed comprising 128 candidate genes (Table 1). For the purposes of variant classification, we considered that genes that had not been previously implicated in the molecular cause of idiopathic or familial OP or OI to be genes of uncertain significance (GUS).

Table 1.

Candidate genes included in the customized targeted gene sequencing panel, and their association to bone fragility

Association with bone fragility	Genes
Associated with nonsyndromic idiopathic or familial OP	DKK1, LRP5, MTHFR, PLS3, WNT1, WNT3A
Associated with Mendelian forms of OP or OI	BMP1, COL1A1, COL1A2, CRTAP, FKBP10, IFITM5, LRP5, NOTCH2, P3H1, P4HB, PLOD2, PPIB, SEC24D, SERPINF1, SERPINH1, SLC34A1, SLC9A3R1, SP7, TMEM38B, WNT1
Associated with other Mendelian diseases with high impact on bone strength, or associated with bone mineral density or fracture risk in major genome-wide association studies (genes of uncertain significance)	ABCF2, ABL1, ADAMTS18, ALDH7A1, ANAPC1, ANOS1, ARHGAP1, AXIN1, C17orf53, C7orf76, CA2, CCDC170, CDC5L, CKAP5, CLCN7, CLDN14, COLEC10, CPED1, CPN1, CREB3L1, CRHR1, CTNNB1, CTSK, CYLD, DCDC1, DCDC5, DHH, DKK1, DLX5, DLX6, DMP1, DNM3, DSPP, ERC1, ESR1, F2, FAM210A, FAM3C, FAM9A, FAM9B, FBN1, FKBP11, FOXC2, FOXL1, FUBP3, GALNT3, GPATCH1, GPR68, HDAC5, IBSP, IDUA, INSIG2, JAG1, KCNMA1, LACTB2, LEKR1, LGR4, LIN7C, LRP4, LRP5, LRP6, MARK3, MEF2C, MEPE, MPP7, NBR1, NTAN1, OSTM1, PDXDC1, PKDCC, PLEKHM1, PLOD3, PTDSS1, PTDSS2, PTHLH, RPS6KA5, RSPO3, RUNX2, SALL1, SHFM1, SLC29A3, SMG6, SMOC1, SNX10, SOST, SOX4, SOX6, SOX9, SP7, SPP1, SPTBN1, STARD3NL, SUCO, SUPT3H, TCIRG1, TMEM263, TNFRSF11A, TNFRSF11B, TNFSF11, TTC21B, USF3, WLS, WNT16, WNT5B, XKR9, ZBTB40, ZNF408

Abbreviations: OI, osteogenesis imperfecta; OP, osteoporosis.

Massively Parallel Sequencing

Targeted regions were captured using a SureSelect kit (Agilent Technologies), following the manufacturer’s instructions. Capture probes were designed using the online Agilent SureDesign software, based on human reference genome GRCh37. All exonic regions and 25-bp boundaries were covered with 3-times minimum tiling. Agilent SureSelect libraries were prepared from 3 μg of genomic DNA sheared using an E220 Focused-Ultrasonicator (Covaris) according to the manufacturer’s instructions. Barcoded libraries were sequenced on a NextSeq500 system (Illumina) and resulting paired-end reads were aligned to the human reference genome GRCh37/hg19 using the Burrows-Wheeler alignment tool. Quality control metrics were obtained through analysis with Qualimap2. Variant calling was performed with Platypus in all the resulting BAM files, and the resulting variants were annotated with ANNOVAR.

Variant filtering was performed to prioritize rare variants affecting protein function, according to these criteria: 1) allele frequency less than 1% both in the Genome Aggregation Database (gnomAD, version 3) and the database of the Online Archive of Brazilian Mutations (ABraOM, Arquivo Brasileiro Online de Mutações, version 1) [35]; and 2) a) variants leading to stop gain, frameshift, in-frame indel, or b) variants predicted to lead to splice site abnormalities according to NNSPLICE, Human Splicing Finder, and NetGene2, or c) nonsynonymous variants predicted to be deleterious according to Sorting Intolerant From Tolerant (SIFT), PolyPhen2, genome evolutionary rate profiling (GERP++), or combined annotation–dependent depletion (CADD).

With regard to the PHRED-like scaled C-scores according to the CADD framework (CADDp), higher scores are associated with variants more likely to be deleterious, but a threshold for deleteriousness is not defined. Scores greater than 10 indicate that the variant is predicted to be among the 10% most deleterious substitutions to the human genome, greater than 20 correspond to the 1% most deleterious, and greater than 30 to the 0.1% most deleterious. In this study, for the purposes of variant prioritization, we have adopted the threshold of greater than 15, as suggested by the developers of CADD and others [36, 37].

Additional data on variants of interest were obtained from ClinVar, the Human Gene Mutation Database, the Leiden Open Variation Database, the Osteogenesis Imperfecta Variation Database, and PubMed. Variants were subsequently classified according to the interpretation criteria proposed by the American College of Medical Genetics and Genomics and the Association for Molecular Pathology (ACMG-AMP) [38].

Copy Number Variation Analysis

Basic copy number variation (CNV) analysis of sequencing data was performed using the COpy Number Targeted Resequencing Analysis (CONTRA) tool. Identified CNVs were verified by BAM file visualization using the Integrative Genomics Viewer (Broad Institute). Verified CNVs were further analyzed according to their size: for regions smaller than 500 bp, a Sanger sequencing-based approach was used, and for larger regions, DNA samples were analyzed by single-nucleotide variation (SNV; formerly SNP array using CytoSNP-850K arrays (Illumina). Briefly, DNA amplification, hybridization, staining, and washing were performed according to the manufacturer’s instructions, and arrays were scanned using the iScan System (Illumina). Raw data were analyzed using BlueFuse Multi, version 1.1 software (Blue Gnome). CNVs confirmed by SNP array were queried in the UCSC genome browser and the Database of Genomic Variants to retrieve previous reports.

Sanger Sequencing

Automated Sanger sequencing was carried out to confirm variants of interest and for segregation analysis when familial samples were available. Polymerase chain reaction–amplified regions (primer sequences are available on request) were purified enzymatically with Illustra ExoProStar (GE Healthcare Life Sciences) and sequenced using a BigDye Terminator, version 3.1 kit (Thermo Fisher Scientific) on an ABI 3130x1 automated DNA sequencer (Thermo Fisher).

Segregation Analysis

First-degree relatives of study individuals were invited to participate. After informed written consent, a detailed medical history was obtained aiming to identify a history of fragility fractures, previous diagnosis of OP/low BMD, or potential causes of secondary OP, and a DXA scan was performed.

Results

Targeted MPS was carried out in 37 participants, including 21 sporadic and 16 familial cases (7 families). For the purposes of genetic analysis, each family was considered as 1 index case, thus rendering a total of 28 index cases. Clinical characteristics of the participants are presented in Table 2, and family pedigrees in Fig. 1. Fifty-four percent of the cohort was male, and the median age at diagnosis was 44 years (range, 19-81 years). Fractures were present in 86% of the cohort, and median age at first fracture was 31.5 years (range, 1-84 years). Fifty-nine percent had at least one vertebral fracture, and 73% one or more nonvertebral fractures.

Table 2.

Clinical characteristics of the 37 individuals analyzed with targeted gene sequencing

ID	Sex	Age at diagnosis, y	BMD Z/T score				Age at first fracture, y	Vertebral fracturesa	Nonvertebral fractures
			Z/T	LS	FN	TH
1	M	57	T	–3.6	–1.2	–0.9	NA	0	0
2	M	77	T	–5.5	–4.1	–3.8	19	3 (M)	1 (ankle)
3	M	39	Z	–4.8	–2.9	NA	6	0	3 (humerus, fibula, scapula)
4	F	33	T	–4.8	–2.8	NA	33	5 (C)	5 (humerus, toe, ribs)
5	M	35	T	–3.5	–2.2	–2.3	8	9 (C)	> 10 (long bones, hand)
6	M	42	Z	–2.7	–2.1	–2.0	7	1 (C)	5 (radius, hand, ribs)
7	F	35	Z	–3.8	–1.4	–1.6	15	0	4 (tibia, hand, foot, shoulder)
8	F	25	Z	–5.7	–4.3	–4.6	7	9 (C)	8 (femur, tibiae, elbows, forearm, clavicle, foot)
9	F	33	Z	–2.3	–2.2	–1.9	33	0b	2 (femur)
10	M	58	T	–4.9	–3.4	–3.5	58	5 (C)	0
11	M	22	Z	–2.1	–0.8	–0.1	22	0	4 (femurs, tibiae)
12	M	44	Z	–4.6	–3.5	–3.5	55	0	2 (mandible, hand)
13	F	38	Z	–3.5	–1.4	–1.6	33	5 (C)	1 (knee)
14	M	64	T	–6.1	–3.8	–4.1	68	11 (M)	5 (femurs, scapula)
15	F	21	Z	–4.4	–4.4	–4.9	1	4 (M)	1 (femur)
16	F	60	T	–4.0	NA	NA	60	4 (C)	5 (humerus, radius, tibia, clavicle, pelvis)
17	F	22	Z	–3.0	–1.5	–1.8	20	1 (M)	5 (radius, ankle, feet, fingers, ribs)
18	M	26	Z	–2.3	–2.6	–2.0	7	0b	3 (radius, tibia, foot)
19	F	49	T	–5.1	–3.1	–3.3	72	4 (M)	1 (forearm)
20	M	45	Z	–3.0	–1.1	–1.5	NA	0	0
21	M	39	Z	–2.8	–0.5	–2.3	39	12 (C)	1 (pelvis)
F1A	M	53	T	–4.5	–4.5	–3.7	2	4 (M)	6 (forearms, ribs)
F1B	F	70	T	–2.7	–2.5	–1.9	60	0	2 (knee, ankle)
F2A	M	52	T	–3.6	–3.0	NA	52	9 (C)	0
F2B	F	49	T	–4.0	–2.1	–1.2	NA	0	0
F3A	F	33	T	–3.7	–2.4	NA	33	4 (C)	0
F3B	F	56	T	–2.4	–1.5	–1.4	56	1 (M)	2 (fibula, ankle)
F4A	F	65	T	–4.2	–3.7	–4.3	30	1 (M)	4 (femur, foot, ribs)
F4B	F	81	T	–1.8	–2.4	–2.7	81	0b	3 (femur, radius, ulna)
F5A	F	65	T	–4.2	–3.8	NA	84	3 (M)	0
F5B	F	54	T	–1.7	–3.0	–3.1	10	0b	2 (forearm, clavicle)
F6A	M	35	Z	–3.5	–1.9	–2.0	19	1 (C)	3 (forearm, humerus, ribs)
F6B	M	32	Z	–2.7	–0.7	–0.2	13	0	3 (tibia, wrist, finger)
F6C	M	61	T	–0.9	–1.3	–0.4	9	1 (M)	4 (forearm, elbow, wrist, foot)
F7A	M	51	T	–3.1	–2.9	–2.8	56	3 (C)	0
F7B	M	19	Z	–2.1	–0.5	–0.9	NA	0	0
F7C	M	20	Z	–3.1	–1.8	–1.6	NA	0	0

Abbreviations: BMD, bone mineral density; F, female; FN, femoral neck; ID, identification; LS, lumbar spine; M, male; NA, not available; TH, total hip.

(C), clinical fractures; (M), morphometric fractures.

No clinical vertebral fractures, but lateral spine radiographs were not available.

Figure 1.

Pedigrees of familial cases. Individuals analyzed by massively parallel sequencing are identified according to their identification in Table 2. In family 4, samples from 2 unaffected family members in generation III were available for segregation analysis; carrier status (c, carrier; nc, noncarrier) are shown for the IDUA p.(His82Gln) and KCNMA1 p.(Arg813Gln) variants, respectively.

The mean coverage of targeted regions ranged from 306 to 1728×, and at least 99.38% of targeted base positions were sequenced at greater than 20×. A total of 1835 variants were called before filtering. After prioritizing rare variants (allele frequency < 1%) with a predicted deleterious impact on the protein, 32 heterozygous variants were identified and subsequently confirmed by Sanger sequencing. Segregation analysis was possible for 19 of these 32 variants and did not support a relevant pathogenic role for 5 variants, which were subsequently excluded. Additional analysis of sequencing data using CONTRA rendered 3 potential CNVs, but only 1 was confirmed by SNP array. This confirmed CNV was also considered a variant of interest, whereas the 2 unconfirmed ones were excluded.

Altogether, 28 allelic variants of interest were identified, but only 3 cases were found to bear pathogenic or likely pathogenic variants (Table 3). Seventeen cases had variants of interest that were classified as variants of uncertain significance (VUS) or were located in GUS (Table 4), mainly because previous reports linking variants or candidate genes to idiopathic or familial osteoporosis are lacking. In 8 cases (Ids 9, 14, 15, 16, 18, F1, F6, and F7), no variants of interest were identified.

Table 3.

Cases with pathogenic or likely pathogenic variants identified by targeted gene sequencing

ID	Variant	Allele frequency		In silico prediction				Previous report [ref.]	ACMG-AMP classification
		gnomADa	ABraOM	SIFTb	PP2c	GERP++d	CADDpe
4	AXIN1 p.(Ala311Thr)	0.0000	0.0000	T	D	5.25	28.4	No	NA (GUS)
	COL1A2 p.(Arg708Gln)	0.0012 (LAT)	0.0016	D	D	5.84	25.1	[40, 48-52]	Pathogenic (II): PS1, PS3, PP3
	FOXC2 p.(Cys498Arg)	0.0034 (ASJ)	0.0008	T	D	4.88	15.5	No	NA (GUS)
	PTDSS2 p.(Val250Met)	0.0060 (ASJ)	0.0041	D	D	4.50	22	No	NA (GUS)
8	NOTCH2 p.(Leu2408His)	0.0032 (NFE)	0.0074	D	D	5.35	25.5	ClinVar	VUS
	PLS3 p.(Leu603Phe)	0.0000	0.0000	D	D	5.59	15.6	No	VUS
	WNT1 p.(Gly169Asp)	0.0009 (AFR)	0.0000	D	D	4.94	35	[39]	Likely pathogenic (III): PS1, PP1, PP2, PP3
F4	IDUA p.(His82Gln)	0.0047 (NFE)	0.0082	T	D	–0.41	9.9	[41, 53, 54]	Likely pathogenic (II): PS3, PM5
	KCNMA1 p.(Arg813Gln)	< 0.0001 (NFE)	0.0000	D	D	4.73	26.9	No	NA (GUS)

Abbreviations: ACMG-AMP, American College of Medical Genetics and Genomics and the Association for Molecular Pathology; CADD, combined annotation–dependent depletion; GERP, genome evolutionary rate profiling; GUS, genes of uncertain significance; ID, identification; SIFT, Sorting Intolerant From Tolerant; VUS, variants of uncertain significance.

Highest ethnicity-specific minor allele frequencies are reported (AFR, African; ASJ, Ashkenazi Jewish; LAT, Latino; NFE, Non-Finnish European).

SIFT output prediction: D, damaging; T, tolerated.

PolyPhen2 HumDiv prediction: B, benign; D, probably damaging; P, possibly damaging.

GERP++ rejected substitutions (RS) score (a suggested threshold for deleteriousness is > 4.4).

PHRED-like scaled C-score according to the CADD framework (a suggested threshold for deleteriousness is > 15); ABraOM, Online Archive of Brazilian Mutations; gnomAD, Genome Aggregation Database; NA, not applicable.

Table 4.

Cases with variants of interest classified as variants of uncertain significance or in genes of uncertain significance

ID	Variant	Allele frequency			In silico prediction			Previous report	ACMG-AMP classification
		gnomADa	ABraOM	SIFTb	PP2c	GERP++d	CADDpe
1	PTDSS1 p.(Gly72Ser)	0.0017 (NFE)	0.0016	D	D	5.54	28	No	NA (GUS)
2	COLEC10 p.(Arg125Trp)	0.0011 (ASJ)	0.0000	D	D	4.1	21.8	No	NA (GUS)
3	COL1A1 p.(Arg528His)	0.0012 (ASJ)	0.0008	D	D	5.06	17.2	ClinVar	VUS
	WLS p.(Gln25His)	0.0011 (LAT)	0.0008	D	D	–2.39	16.3	No	NA (GUS)
	WNT1 p.(Cys93Tyr)	0.0000	0.0000	D	D	5.04	24.7	No	VUS
5	P3H1 p.(Pro358Thr)	< 0.0001 (LAT)	0.0000	D	P	5.84	15.2	No	VUS
6	GPR68 exonic deletion	NA	NA	NA	NA	NA	NA	No	NA (GUS)
7	NBR1 p.(Gly759Val)	0.0001 (NFE)	0.0000	D	D	5.64	25.8	No	NA (GUS)
10	CCDC170 p.(Glu451Lys)	0.0025 (AFR)	0.0000	D	D	1.21	33	No	NA (GUS)
11	LACTB2 p.(144_147del)f	0.0051 (NFE)	0.0057	NA	NA	NA	NA	No	NA (GUS)
	NOTCH2 p.(Leu2408His)f	0.0032 (NFE)	0.0074	D	D	5.35	25.5	ClinVar	VUS
12	PKDCC p.(Asn210Lys)	0.0000	0.0000	D	D	3.12	33	No	NA (GUS)
13	PLS3 p.(Arg94Cys)	0.0000	0.0020	D	D	4.00	34	No	VUS
17	AXIN1 p.(Val683Met)	< 0.0001 (SAS)	0.0000	T	D	2.18	22.7	No	NA (GUS)
19	NBR1 p.(Asp40Gly)	0.0000	0.0000	T	D	5.26	23.5	No	NA (GUS)
20	USF3 p.(Ser1425Leu)	0.0005 (LAT)	0.0000	D	B	5.06	22.3	No	NA (GUS)
21	LACTB2 p.(144_147del)f	0.0051 (NFE)	0.0057	NA	NA	NA	NA	No	NA (GUS)
F2	GALNT3 p.(Lys347Thr)	0.0000	0.0000	T	B	4.41	18.5	No	NA (GUS)
	WNT1 p.(Thr336Met)	0.0093 (AFR)	0.0017	D	D	3.5	19.6	No	VUS
F3	BMP1 p.(Arg371His)	0.0062 (NFE)	0.0066	T	D	5.77	26.4	ClinVar	VUS
F5	ANAPC1 p.(Val1052Met)	0.0001 (LAT)	0.0000	D	B	4.58	24.1	No	NA (GUS)

Abbreviations: ACMG-AMP, American College of Medical Genetics and Genomics and the Association for Molecular Pathology; CADD, combined annotation–dependent depletion; GERP, genome evolutionary rate profiling; GUS, genes of uncertain significance; ID, identification; SIFT, Sorting Intolerant From Tolerant; VUS, variants of uncertain significance.

Highest ethnicity-specific minor allele frequencies are reported (AFR, African; ASJ, Ashkenazi Jewish; LAT, Latino; NFE, Non-Finnish European; SAS, South Asian).

SIFT output prediction: D, damaging; T, tolerated.

PolyPhen2 HumDiv prediction: B, benign; D, probably damaging; P, possibly damaging.

GERP++ rejected substitutions (RS) score (a suggested threshold for deleteriousness is > 4.4).

PHRED-like scaled C-score according to the CADD framework (a suggested threshold for deleteriousness is > 15); ABraOM, Online Archive of Brazilian Mutations; gnomAD, Genome Aggregation Database; NA, not applicable.

Recurring variants.

One pathogenic variant, COL1A2 p.(Arg708Gln), was identified in patient 4 in combination with 3 VUS (Table 3). Patient 4 is a woman who, at the time of menopause (early, at age 33 years), already had multiple vertebral and nonvertebral fractures and low bone mass (see Table 2). Two variants were classified as likely pathogenic: WNT1 p.(Gly169Asp), found in association with 2 VUS in patient 8, and IDUA p.(His82Gln), found in association with 1 VUS in family 4.

In particular, patient 8, a woman with multiple long bone fractures since childhood and severe vertebral fractures during labor but no extraskeletal features suggestive of OI, carried an association of 3 heterozygous variants: WNT1 p.(Gly169Asp), PLS3 p.(Leu603Phe), and NOTCH2 p.(Leu2408His) (Fig. 2). The likely pathogenic WNT1 p.(Gly169Asp) variant has been reported in 2 Chinese individuals with OI in compound heterozygosity with other WNT1 variants [39], but patient 8 does not bear additional WNT1 variants. The remainder variants are currently considered VUS. Segregation analysis has shown that the mother of patient 8 carries the WNT1 and PLS3 variants. She was diagnosed with postmenopausal OP with an LS and FN T score of –4.4 at age 62 years and has a history of multiple fragility fractures of tibiae and vertebrae since age 45 years. The father, carrier of the NOTCH2 variant, was also diagnosed with OP at age 67 years, with an FN T score of –3.1 and no clinical fractures. The association of these 3 variants could be seen to confer the more severe phenotype of patient 8.

Figure 2.

Patient 8’s clinical and molecular characterization. A, Magnetic resonance imaging of patient 8, a 25-year-old woman with early-onset osteoporosis (OP), shows several thoracic and lumbar vertebral fractures. B, Her family pedigree shows that her mother also has OP and fractures, and her father has OP without fractures. Other cases of bone fragility were identified in her mother’s family. Dual-energy x-ray absorptiometry (DXA) scans (table) show variable degrees of bone loss in lumbar spine (LS), total hip (TH), and femoral neck (FN) in patient 8, and her mother and father. C, Electropherograms of the 3 heterozygous variants identified in patient 8; her mother carries the PLS3 and WNT1 variants, whereas her father carries the NOTCH2 variant.

Altogether, variants of interest were found in 23 candidate genes (18% of gene panel). Seven of these genes had previous association with idiopathic or Mendelian forms of bone fragility (BMP1, COL1A1, COL1A2, NOTCH2, P3H1, PLS3, and WNT1) and the remaining 16 were GUS (ANAPC1, AXIN1, CCDC170, COLEC10, FOXC2, GALNT3, GPR68, IDUA, KCNMA1, LACTB2, NBR1, PKDCC, PTDSS1, PTDSS2, USF3, and WLS). Six variants have been previously reported in ClinVar or in the literature [39-41] (Tables 3 and 4).

Among the GUS bearing variants of interest, NBR1 and GPR68 are noteworthy. Variants in NBR1 were identified in 2 cases: patient 7, a 35-year-old woman with an LS Z score of –3.8 and multiple fractures, bearing the NBR1 p.(Gly759Val) variant, and patient 19, a 49-year-old woman with an LS T score of –5.1 shortly after menopause, carrying the NBR1 p.(Asp40Gly) variant. The p.(Asp40Gly) variant was also found in patient 19’s 50-year-old son, who has an LS T score of –3.9, but not in her 42-year-old son with an LS Z score of –0.6, corroborating a potential association with bone fragility. As for GPR68, a previously unreported heterozygous deletion of its single coding exon was found in patient 6, an otherwise healthy man diagnosed with OP at age 42 years after a vertebral fracture during light physical activity (see Table 2). Subsequent analysis with SNP array confirmed this CNV, which spanned 41 543 bp comprising exon 2 of GPR68 and the last 3 exons of neighboring gene DGLUCY.

Discussion

In this study, we aimed to analyze the genetic contribution to unusual cases of OP. From a bone fragility perspective, these cases represent a phenotype of intermediate severity between common postmenopausal/senile OP and rare early-onset Mendelian disorders with multiple fractures, such as OI. Individuals in our cohort had a younger age at diagnosis and higher fracture prevalence than what would typically be expected for common OP, whereas syndromic or extraskeletal features suggesting Mendelian disorders were absent. Therefore, hoping to broaden our knowledge on the genetic architecture of OP, a targeted sequencing panel was customized to include a large number of candidate genes for bone fragility to be queried in this cohort. Although we have identified 28 rare copy-number or nonsynonymous variants predicted to alter protein function in 71.4% of the cohort, only 3 were classified as pathogenic or likely pathogenic.

Attributing pathogenicity to genetic variants identified by MPS is challenging, particularly when studying a common disorder such as OP. For Mendelian disorders, attributes such as absence in population databases, high impact on protein structure, and previous knowledge of loss of function as a disease mechanism can help determine pathogenicity. The paucity of studies investigating idiopathic, severe, or familial OP means that gene loss of function will not be clearly defined as a disease mechanism for most candidate genes, and that previous descriptions of identified variants will generally be lacking. As a means of identifying potentially pathogenic variants, we have prioritized rare variants with predicted deleteriousness. In silico tools are widely used for these purposes and include function prediction scores (eg, SIFT and PolyPhen2), conservation scores (eg, GERP++), and integrative annotation scores (eg, CADD). Though useful for prioritization, in silico prediction evidence serves only as support for pathogenicity, rather than determination [37, 38, 42].

Although almost all variants identified in this study were nonsynonymous substitutions, a large deletion of the single coding exon of GPR68 was identified in patient 6. Candidate gene GPR68 was included in the sequencing panel because of its proximity to the strong GWAS locus rs1286083 [22]. In vitro studies have supported a role for GPR68 in osteoclastogenesis, and Gpr68-deficient mice have reduced osteoclast differentiation [43, 44]. Recently, homozygous defects in this gene have been associated with human amelogenesis imperfecta, a disease with altered enamel mineralization [45].

Twenty-five percent of variants were identified in genes previously associated with severe OP or to mild OI, such as WNT1, PLS3, COL1A1, and COL1A2. Notably, despite high sequencing coverage only one variant of interest was found in LRP5, which was subsequently excluded by segregation analysis. Even though this gene is frequently associated with idiopathic OP [32, 40], a low prevalence of LRP5 variants has also been reported in a Belgian cohort of idiopathic OP in men and a British cohort of juvenile idiopathic OP [31, 46].

Interestingly, 2 distinct variants were found in the GWAS candidate gene NBR1 in patients 7 and 19. Transgenic Nbr1 mouse models have been reported with altered BMD, abnormal osteoclast physiology, and osteoblast differentiation, as well as changes in bone microarchitecture [47]. Collectively, these data warrant further exploration of NBR1 involvement with bone strength in humans.

Our study has limitations. GWAS identify loci associated with low BMD or fractures, and candidate genes derived from these studies were based on criteria such as proximity to strongest signaling SNV and available data on gene function and animal models; however, they may not necessarily represent the biologically relevant candidate genes in relation to those loci. The selected cohort was heterogeneous, with a wide range of age at diagnosis and age at first fracture; nevertheless, the individuals are representative of the phenotype we aimed to investigate. Even though a broad selection of candidate genes was analyzed, we could not find variants of interest in 28.6% of the cohort. In particular, no variants were identified in 3 out of the 7 families studied (families 1, 6, and 7), which had affected individuals in multiple consecutive generations, yielding a high suspicion of an underlying genetic cause. Considering the high quality and depth of our gene panel sequencing, most likely the gene or genes associated with the disease in these families were not included in our panel, and, therefore, exome or genome sequencing might be able to determine their genetic etiology in the future.

As discussed earlier, attribution of pathogenicity to variants identified by MPS is not straightforward. In this regard, applying the ACMG-AMP proposed criteria for interpretation of sequence variants to the 28 variants, we have identified yielded 25 VUS, 2 likely pathogenic variants, WNT1 p.(Gly169Asp) and IDUA p.(His82Gln), and 1 pathogenic variant, COL1A2 p.(Arg708Gln). The IDUA p.(His82Gln) and COL1A2 p.(Arg708Gln) variants have been previously reported in ClinVar, with an interpretation of pathogenicity discrepant to ours, given that they had been studied in different conditions (mucopolysaccharidosis type I and OI, respectively). Nevertheless, and corroborating our interpretation of pathogenicity of these variants, COL1A2 p.(Arg708Gln) has been previously associated with severe OP, and for both variants a functional impact on encoded proteins has been demonstrated in the literature [40, 41, 48]. The majority of GUS in our panel and the scarcity of published genetic analyses of idiopathic, familial, or severe OP explain the high proportion of VUS. Although the ACMG-AMP guidelines are invaluable in the diagnosis of classic Mendelian disorders in clinical practice, they are not primarily intended for the analysis of variants in multigenic non-Mendelian complex disorders such as osteoporosis. Nevertheless, it is telling that in a cohort enriched for extreme presentations of OP only a few pathogenic genetic variants could be identified. Future endeavors for the in vitro functional characterization of likely pathogenic variants and VUS identified in this cohort would be of great value to gauge their pathogenicity.

In conclusion, targeted MPS of 128 candidate genes in 28 cases of idiopathic, severe, or familial OP rendered the identification of pathogenic or likely pathogenic variants in only 3 cases (11%). Established candidate genes such as WNT1, PLS3, COL1A1, and COL1A2 were validated in association with OP in our cohort, and a potential additional line of evidence has arisen for GWAS candidate genes GPR68 and NBR1, warranting further studies. Although we hope that continuing efforts to identify genetic predisposition to OP will lead to improved and personalized care in the future, it is important to recognize that the likelihood of identifying actionable pathogenic variants in intriguing cases of idiopathic or familial osteoporosis is currently low.