
==== Front
J Bone Miner Res
J Bone Miner Res
jbmr
Journal of Bone and Mineral Research
0884-0431
1523-4681
Oxford University Press

39126371
10.1093/jbmr/zjae122
zjae122
Research Article
AcademicSubjects/MED00010
AcademicSubjects/MED00160
AcademicSubjects/MED00250
Gene-based association analysis of a large patient cohort provides insights into genetics of atypical femur fractures
https://orcid.org/0000-0002-6748-6303
Zhou Wei Department of Internal Medicine, Erasmus MC, University Medical Center Rotterdam, 3015 GD Rotterdam, The Netherlands

Ås Joel Department of Medical Sciences, Uppsala University Hospital, Uppsala 75185, Sweden

Shore-Lorenti Catherine Department of Medicine, School of Clinical Sciences, Monash University, Clayton VIC 3168, Australia

https://orcid.org/0000-0002-8846-6168
Nguyen Hanh H Department of Medicine, School of Clinical Sciences, Monash University, Clayton VIC 3168, Australia

van de Laarschot Denise M Department of Internal Medicine, Erasmus MC, University Medical Center Rotterdam, 3015 GD Rotterdam, The Netherlands

Sztal-Mazer Shoshana Department of Endocrinology and Diabetes, The Alfred Hospital, Melbourne VIC 3004, Australia
Department of Public Health and Preventative Medicine, Monash University, Melbourne VIC 3004, Australia

Grill Vivian Department of Endocrinology and Diabetes, Western Health, Melbourne VIC 3011, Australia

Girgis Christian M Department of Diabetes and Endocrinology, Westmead Hospital, Westmead, NSW 2145, Australia
Faculty of Medicine and Health, The Sydney University, Camperdown NSW 2050, Australia

Stricker Bruno H Ch Department of Epidemiology, Erasmus MC, University Medical Center Rotterdam, Rotterdam 3015 GD, The Netherlands

van der Eerden Bram C J Department of Internal Medicine, Erasmus MC, University Medical Center Rotterdam, 3015 GD Rotterdam, The Netherlands

Thakker Rajesh V Academic Endocrine Unit, Radcliffe Department of Medicine, University of Oxford, Oxford OX3 9DU, United Kingdom
National Institute for Health Research Oxford Biomedical Research Centre, Oxford OX3 9DU, United Kingdom

https://orcid.org/0000-0001-5035-127X
Appelman-Dijkstra Natasha M Department of Internal Medicine, division endocrinology, Leiden University Medical Center, Leiden 2333 ZA, The Netherlands

Wadelius Mia Department of Medical Sciences, Uppsala University Hospital, Uppsala 75185, Sweden

Clifton-Bligh Roderick J Faculty of Medicine and Health, The Sydney University, Camperdown NSW 2050, Australia

Hallberg Pär Department of Medical Sciences, Uppsala University Hospital, Uppsala 75185, Sweden

Verkerk Annemieke J M H Department of Internal Medicine, Erasmus MC, University Medical Center Rotterdam, 3015 GD Rotterdam, The Netherlands

van Rooij Jeroen G J Department of Internal Medicine, Erasmus MC, University Medical Center Rotterdam, 3015 GD Rotterdam, The Netherlands

Ebeling Peter R Department of Medicine, School of Clinical Sciences, Monash University, Clayton VIC 3168, Australia

https://orcid.org/0000-0001-9186-3423
Zillikens M Carola Department of Internal Medicine, Erasmus MC, University Medical Center Rotterdam, 3015 GD Rotterdam, The Netherlands

Corresponding author: M. Carola Zillikens, Department of Internal Medicine, Erasmus MC, University Medical Center Rotterdam, Dr. Molewaterplein 40, 3015 GD Rotterdam, The Netherlands (m.c.zillikens@erasmusmc.nl).
9 2024
10 8 2024
10 8 2024
39 9 13151326
18 12 2023
18 6 2024
08 7 2024
23 8 2024
© The Author(s) 2024. Published by Oxford University Press on behalf of the American Society for Bone and Mineral Research.
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (https://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com

Abstract

Several small genetic association studies have been conducted for atypical femur fracture (AFF) without replication of results. We assessed previously implicated and novel genes associated with AFFs in a larger set of unrelated AFF cases using whole exome sequencing (WES). We performed gene-based association analysis on 139 European AFF cases and 196 controls matched for bisphosphonate use. We tested all rare, protein-altering variants using both candidate gene and hypothesis-free approaches. In the latter, genes suggestively associated with AFFs (uncorrected p-values <.01) were investigated in a Swedish whole-genome sequencing replication study and assessed in 46 non-European cases. In the candidate gene analysis, PLOD2 showed a suggestive signal. The hypothesis-free approach revealed 10 tentative associations, with XRN2, SORD, and PLOD2 being the most likely candidates for AFF. XRN2 and PLOD2 showed consistent direction of effect estimates in the replication analysis, albeit not statistically significant. Three SNPs associated with SORD expression according to the GTEx portal were in linkage disequilibrium (R2 ≥ 0.2) with an SNP previously reported in a genome-wide association study of AFF. The prevalence of carriers of variants for both PLOD2 and SORD was higher in Asian versus European cases. While we did not identify genes enriched for damaging variants, we found suggestive evidence of a role for XRN2, PLOD2, and SORD, which requires further investigation. Our findings indicate that genetic factors responsible for AFFs are not widely shared among AFF cases. The study provides a stepping-stone for future larger genetic studies of AFF.

whole exome sequencing
atypical femur fractures
bisphosphonates
osteoporosis
gene
Jaap Schouten Foundation 10.13039/501100018722 National Health and Medical Research Council of Australia GNT1197958 Wellcome Trust Investigator Award 106995/Z/15/Z National Institute for Health Research 10.13039/501100000272 NIHR Senior Investigator Award NF- SI-0514–10091 Netherlands Consortium for Healthy Aging 050-060-810
==== Body
pmcIntroduction

An atypical femur fracture (AFF) is considered a devastating rare adverse effect of antiresorptive therapy, such as bisphosphonates (BP), next to osteonecrosis of the jaw. Despite being highly effective in preventing osteoporotic fractures, BP treatment has decreased due to concerns among patients and physicians about these adverse effects.1,2 Therefore, a better understanding of AFF is needed to help recognize patients at risk and guide the use of antiresorptive therapy.

AFF has an incidence of 3-10 per 100 000 person-years in the general population with a duration-dependent association between BP use and AFF risk.3 Patients with more than 8 yr of BP use have 43-fold increased risk compared with those with less than 3 mo of use.4 AFFs have distinct radiological features, as defined in the Task Force report of the American Society for Bone and Mineral Research (ASBMR), in comparison to plain osteoporotic fragility fractures.5 Although rare, these fractures have a significant impact on patients due to an increased risk of delayed healing and high chance of contralateral femur fractures, and reduced mobility and level of function.6

So far, it is not understood why only a small fraction of patients with prolonged BP use develops AFF. One possible explanation is a genetic susceptibility to AFF, supported by the occurrence of multiple AFF cases in several families.7-10 Moreover, our systematic review of AFF cases with monogenic bone disorders suggests that genetic factors related to these disorders may predispose to AFF.8 In line with this, we showed that 15% of Dutch AFF patients have pathogenic variants or likely pathogenic variants in genes related to monogenic bone disorders.11 In other studies, it was suggested that genes related to abnormal bone mineralization are involved in the predisposition to AFF.12,13 Additionally, genes related to the mechanisms of BP action have been identified as potential contributors to AFF. For example, variants in GGPS1 and CYP1A1 genes have been associated with AFF, possibly via altering BP action in carriers of these variants.9,14,15 Another study discovered a variant in the ATRAID gene in some AFF patients, which influenced cellular sensitivity to BPs.16

To identify genetic variants associated with non-familial AFF, two exon-wide and genome-wide single variant association analyses have been performed using a small number of AFF cases.17,18 However, neither study has established an association between a specific genetic variant and AFF, and none of the reported findings have been replicated to date.

In the current study, we conducted gene-based association analyses focusing on rare, protein-altering variants. We used a gene-based approach instead of a single-variant based approach because the latter requires a larger sample size and number of carriers per variant to demonstrate a statistically significant effect.19 The first aim of the study was to assess the association between AFF and candidate genes comprising: (1) genes responsible for monogenic bone disorders including those involved in bone mineralization; (2) genes involved in the mechanism of BP action; and (3) genes implicated in previous AFF genetic studies. The second aim of the study was to explore novel genes for AFF using a hypothesis-free approach.

Materials and methods

Discovery dataset

In total, 185 adult multi-ethnic, unrelated AFF patients from the Netherlands (n = 74), Singapore (n = 2), and Australia (n = 109) were included in the study between September 2013 and January 2021. Among the included AFF patients, 91% had been exposed to BPs before their AFF. All patients signed written informed consent to participate in this study. The phenotypic data were extracted from medical files and/or from questionnaires. The radiological features of all the included AFF patients fulfilled the revised case definition for AFF published in the second ASBMR Task Force Report in 2014.5 The study was approved by the Medical Ethical Committee (METC) of Erasmus MC under number MEC-2013-264 and Melbourne Health Human Research Ethics Committee (HREC) (HREC/14/MH/160) in Australia.

Controls were from the Rotterdam Study (RS), which is a prospective cohort study of persons living in the well-defined Ommoord district in the city of Rotterdam in the Netherlands who were ≥55 yr old at the start of the study in 1990.20 The controls included in this study were part of the original cohort (RS1) recruited during 1989 and 1993 who were randomly selected for whole exome sequencing.21 After WES data quality control, 2604 samples remained in this dataset.21 Ethics statement for this study is included in Supplementary methods. The diagnosis of osteopenia and osteoporosis was obtained from baseline interviews or defined using measurements of dual-energy X-ray absorptiometry scanning at their follow-up visits (Bone mineral density T-score ≤ −1 for osteopenia and ≤−2.5 for osteoporosis). The medication history was extracted from pharmacy records using anatomical therapeutical chemical code. Their history of AFF was not available, and we assumed that none of the controls had experienced an AFF due to its rarity in the general population.

WES and data processing

WES was performed on blood-derived DNA samples using Illumina paired-end sequencing. Cases were sequenced using the EZ MedExome Capture kit (Roche Nimblegen, Inc, Madison, WI, USA) and control samples were sequenced using the Nimblegen SeqCap EZ V2 kit (Roche Nimblegen, Inc., Madison, WI, USA). The raw sequencing data were pre-processed using a standard pipeline as previously described11 and aligned to the human reference genome UCSC build hg19. Details of quality control are included in Supplementary methods.

Because different exome capture kits were used to sequence case and control samples, a stringent quality control at the variant level was applied to harmonize the two datasets. After combining genomic variant call format (VCF) files of both case and control samples in one VCF file, genotypes with DP (total depth) < 10, or GQ (genotype quality) < 20 were marked as missing. Following that, we removed variants with missing rates >90% across all samples. Furthermore, we compared the proportions of missing genotypes per variant between case samples and control samples using binomial tests as described by Raghavan and colleagues previously.22 With this method, we removed variants with differential genotype missing rates between the two groups (p-value <.05 in a binomial test). Finally, variants with QD (quality of depth) score <5 were removed.

Variant annotation

Variants were annotated using ANNOVAR (version 2019-10-24). Variant frequencies were obtained from gnomAD v2.1.23 Combined Annotation Dependent Depletion (CADD) scores downloaded from ANNOVAR were used to predict the likelihood of a variant having a deleterious impact on protein function.24 Single nucleotide variant (SNV) CADD scores of 10 and 20 indicate that the variants are among the top 10% and 1% damaging SNVs, respectively. All details of bioinformatic pre-processing are described in Supplementary methods.

Gene-based association analyses

The gene-based association study was performed using both a candidate gene approach and a hypothesis-free approach. The latter was followed by downstream analyses and replication in another dataset for prioritization of genes with suggestive signals. The flowchart of the analyses is shown in Figure 1.

Figure 1 The flowchart of the gene-based association analyses using a candidate gene approach and a hypothesis free approach. Abbreviations: OR, odds ratio. CADD, combined annotation dependent depletion score. SKAT-O test, statistical method that optimally combines kernel and burden tests.

Samples included for association analyses

Principal component analysis (PCA) was performed using WES data to determine European and Asian origin as described in Supplementary methods.25 In total, 139 of the 185 included unrelated AFF cases and 2588 of the 2604 controls were of European origin and were included for association analysis. Non-European cases were used for comparison of carrier frequencies of the genes with a suggestive signal, because no controls were available for these ethnicities. The first analyses involved comparing 139 European AFF cases with 196 European controls matched for BP use. Since most AFF patients had prolonged BP treatment, the matched controls were arbitrarily defined as controls that had been exposed to BPs for at least 3 yr. The percentages of females, individuals with a history of osteoporosis or osteopenia, and individuals with GC use over 1 yr were similar between 139 European AFF cases and 196 BP-matched controls (Table 1).

Table 1 Characteristics of 139 AFF cases and 196 BP-matched controls of European origin.

	AFF patients
(n = 139)	BP-matched controls
(n = 196)	p-value	
Age of AFF onset, yr (Mean, SD)	68 ± 13	N/A	N/A	
Female (N, %)	121 (87.1%)	167 (85.2%)	>.05	
BP use (N, %)	126 (90.6%)	196 (100.0%)	<.05	
BP duration, yr (Mean, SD)	8.7 ± 5.7	6.9 ± 3.0	<.05	
GC use > 1 yr (Mean, SD)	56 (40.3%)	91 (46.4%)	>.05	
OP history (N, %)	116 (83.5%)	178 (90.8%)	>.05	
Abbreviations: AFF, atypical femur fractures; BP, bisphosphonate; GC, glucocorticoid use for over one year; OP, osteopenia or osteoporosis.

Variant filtering from WES data

For both the candidate gene and the hypothesis-free analysis, we included rare, protein-altering variants. Variants were considered rare if they had an allele frequency < 1% in the overall population of the gnomAD database. The protein-altering variants include nonsynonymous, stopgain, stoploss variants, exonic insertions and deletions (indels), and variants residing in splicing regions.

Statistical models

The associations were evaluated using the SKAT-O test, which optimally combines kernel and burden tests to maximize power.26 The SKAT-O test was performed using the “SKATbinary” function implemented in the R package SKAT (version 2.2.4).27 Genotypes were coded additively as 0, 1, or 2. For SKAT-O analysis on X chromosome data from males, homozygous and hemizygous genotypes were coded as 0 and 2, respectively. The associations were first assessed without including covariates and genes. Associations with nominal p-values <.01 were further characterized by correcting for sex and the first four principal components (PCs).28 PCA was performed again on the WES data, but this time only using the sample of European origin, and the first four largest PCs were used to adjust for population stratification within the European dataset. A suggestive signal was defined as having nominal p-values <.01 in both tests. False discovery rate (FDR) was calculated by the Benjamini–Hochberg procedure to control for type-I error rate. FDR < 0.05 was considered statistically significant. Because the SKAT-O tests do not provide effect estimates, we also used the Combined Multivariate and Collapsing method to estimate odds ratio (OR) using a two-sided Fisher’s exact test.19 The estimates were obtained from the “fisher.test” function provided in R package stats (version 4.0.0). Genes with at least one case or control carrying rare, protein-altering variant(s) were included for analyses. This approach compares the number of carriers between cases and controls in a certain gene.

Candidate gene analysis

An extensive list of genes was selected for the candidate gene analysis, as presented in Table S1. The list includes genes associated with monogenic bone disorders or bone mineralization and BP metabolism, and genes reported in previous family-based genetic studies and genetic association studies of AFF. The selection of candidate genes is described in detail in Supplementary methods. For each gene, first, gene-based association analysis was conducted on all the rare (allele frequency < 1%), protein-altering variants (CADD ≥0). Additionally, genes were analysed separately on rare, protein-altering variants with a CADD score ≥ 10 and those with a CADD score ≥ 20. If pathogenic variants in a gene are associated with AFF, we expect to observe higher ORs when restricting the analysis to increasingly damaging variants (CADD ≥ 0, CADD ≥ 10, CADD ≥ 20). Therefore, the criteria for a suggestive signal associated with an increased risk of AFF were defined as: (1) OR ≥ 1 at CADD ≥ 0; (2) an increase in ORs when restricting the analysis to increasingly damaging variants (CADD ≥ 0, CADD ≥ 10, CADD ≥ 20); and (3) a p-value <.01 in any results using SKAT-O.

Hypothesis-free analysis

Under a hypothesis-free approach, association analyses were performed for all the genes that have at least one case or control carrying rare (allele frequency < 1%), protein-altering variant(s) after WES variant filtering.

Additional quality check

The presence and quality of the indels in the suggestive signals were manually reviewed using Integrative Genomics Viewer (IGV version 2.11.9), which displays the alignment of sequence reads in the region (from binary alignment map files). Indels regarded false positive in IGV are those located in low-quality regions (low coverage, excessive low-quality reads, strand bias, highly polymorphic, etc.).29,30 Signals driven by false positive indels were removed from downstream analyses.

Downstream analysis prioritizing the 10 suggested genes

Analysis using unmatched controls

To obtain more robust results, we increased sample size by including unmatched control samples for analysis. The genes with suggestive signals in the hypothesis-free analysis (p-value <.01) were tested for their association with AFF by comparing 139 European cases with all 2588 European controls regardless of BP use. The results of the analysis comparing cases with controls unmatched for BP use (hereafter referred to as unmatched controls) may be biased due to differences in the prevalence of the indications of BP use between the two groups. Despite this, if a gene still appeared to be significant when comparing cases with unmatched controls, it was considered more likely to be a true positive finding.

Prevalence of rare variants in non-European AFF cases

Gene-based association analysis could not be conducted for non-European AFF cases due to the unavailability of non-European controls. We present the percentage of carriers of rare, protein-altering variants in the suggested genes in Asian cases and cases of other ethnicities. A similar or higher prevalence of cases with rare, protein-altering variants in a gene compared with that in the European cases was considered suggestive of a true association with AFF.

Comparison of variant CADD scores between cases and controls

Rare, protein-altering variants relevant to AFF were expected to have a higher CADD score compared with variants with small impact (lower CADD scores) on the same gene. To make a comparison between the CADD scores of the variants carried by cases and controls when CADD score annotation was missing for a variant, the median scores for the variant type were used. This resulted in a score of 37 for frameshift indels and stop-gain variants and 15 for in-frame indels, as reported previously.31 Genes showing higher CADD scores for variants in cases than in controls are considered more likely to be true positive results.

Comparison with previous genetic studies of AFF

The candidate gene analysis described above did not include the 132 genes reported in Table 2 of the paper by Garcia-Giralt et al.32 due to its extensive length and lack of controls for association analysis. However, the genes with suggestive signals in the hypothesis-free analysis were compared with these 132 genes. We also investigated whether these genes were associated with the suggestive signals (p-value <1e−5) reported in the GWAS by Kharazmi et al.18 or their LD blocks (R2 > 0.2) by utilizing expression quantitative trait loci (eQTL) results from the GTEx Portal (Release V8, https://gtexportal.org/).

Table 2 Number of variants and genes remaining in three groups, filtered by increasingly stringent criteria.

Steps	Three groups of variants	Number of variants	Number of genes	
1	All rare, protein-altering variants	511	148	
2	Rare, protein-altering variants with CADD ≥10	407	139	
3	Rare, protein-altering variants with CADD ≥20	186	96	

Gene annotation

The suggested genes were annotated with multiple sources to identify bone-related functions, including (1) genome-wide association studies (GWAS), (2) International Mouse Phenotyping Consortium (IMPC) data, (3) Human Phenotype Ontology, (4) Gene ontology (GO) biological process pathways, (5) the ExAC gene constraint scores, (6) the osteoblasts and osteoclasts expression database, and (7) the bone biopsy expression database (Osteogene).33 Detailed description of these datasets can be found in the Supplementary methods.

Replication dataset

The genes identified by the discovery data were analyzed in individuals from the Swedegene project, which is a Swedish national biobank of adverse drug reactions.34 Whole-genome sequencing data were available for 834 individuals. A total of 53 cases (50 female, 3 male) had AFF associated with BP. The remaining 781 controls had experienced different adverse reactions to a multitude of drugs. Clinical data (Drugs treatment history, demographics, laboratory data, and ancestry) were obtained through interviews, medical records, and questionnaires. All 53 cases of AFF fulfilled the revised case definition for AFF published in the second ASBMR Task Force Report in 2014. The study was approved by the regional ethical review board in Uppsala (2010/231 Uppsala).

In the questionnaire the parents’ place of birth was used as a proxy for genetic ancestry. Among the cases 49 have Swedish origin, 2 from the Nordic countries, 1 from Oceania, and 1 with missing data. Among the unmatched controls, 87% had both parents born in Sweden and 96% had at least one parent with Sweden as birth place. Of the remaining, 20 individuals were missing data representing 3%. There were 39 controls matched for a current diagnosis of osteoporosis and/or BP/denosumab use (n = 39) without a current cancer diagnosis. However, the treatment time is unknown. Among these, 38 are of Swedish origin and 1 is of Nordic origin.

DNA extraction was performed using peripheral blood and DNA libraries were prepared using TrueSeq PCRfree preperation kits, with a target size of 350 bp. Pair-end sequencing was performed on Illumina HiSeq X (sequencing chemistry 2.5) with a read-length of 150 bp to an average coverage of 30. Samples were aligned using BWA-MEM 0.7.1235 to Genome Reference Consortium Human Build 37 (GRCh37) and index using samtools 0.1.19.36 Local realignment was performed using Genome Analysis Toolkit (GATK) 3.3 RealignerTargetCreator and GATK 3.3 IndelRealigner.37 Picard MarkDuplicates 1.120 were used for PCR deduplication and base quality recalibration tables were generated using GATK 3.3 BaseRecalibrator. These steps were performed by the Swedish National Genomics Infrastructure in Uppsala.

Following the GATK best practice workflow,38 per sample variant calling was performed using GATK 3.8 Haplotypecaller39 and joint genotyping was performed on multi-sample 10 Mb subsections of all patients within the cohort using GATK 3.8 GenotypeGVCF. The indel and SNP variants scores were recalibrated using GATK Variant Quality Score Recalibration.

The subsequent replication dataset was analyzed with the same QC steps and statistical methods. Cases were compared with the 781 unmatched controls as well as to the subset of 39 matched controls.

Results

Candidate gene analysis

We performed a candidate gene analysis on 139 European AFF cases and 196 matched European controls, selecting a total of 203 genes for analysis (Table S1). Variants in these genes were filtered and analysed in three variant groups: (1) all rare, protein-altering variants; (2) rare, protein-altering variants with CADD ≥10; and (3) rare, protein-altering variants with CADD ≥20. Across these three groups, an increasingly stringent criterion was applied, with fewer genes remained in the analysis (Table 2). The SKAT-O and Fisher’s exact tests results are presented in Table S2. By SKAT-O correcting for sex and the four largest PCs, PLOD2 was the only gene showing a p-value below .01 in the analyses of variants with CADD ≥ 0 and CADD ≥ 10, but not in the analysis of variants with CADD ≥ 20. None of the other genes had a p-value below .01. However, PLOD2 showed similar effect estimates in all three analyses with nested CADD scores, where an increasing estimated effect might be expected if variants with a higher CADD score are more likely to be damaging (and thus more strongly overrepresented in cases over controls). Eleven other candidate genes had a trend of increasing ORs across the three variant groups (OR > 1 in the first group, highest OR in the last group) (Figure S1). In comparison, a similar number of genes had a trend of decreasing ORs across the three variant groups (OR < 1 in the first group, lowest OR in the last group) (Figure S2). This suggests that both of these results are likely to be chance findings, especially since none were statistically or suggestively significant.

Hypothesis-free analysis

For the hypothesis-free analysis, 139 European cases and 196 matched European controls were compared. After filtering of WES data, 43 101 variants in 13 198 genes were included for gene-based association analyses using SKAT-O. Without correcting for covariates, genetic inflation of the analysis model was λ = 1.08 (QQ-plot: Figure S3). In total, 57 genes had nominal p-values less than .01 (Table S2), but none of the associations were statistically significant after correcting for multiple testing (data not shown). After correcting for sex and the four largest PCs, 14 genes showed suggestive signals (p-value <.01). Four of these gene signals were driven by an indel that was regarded as likely false positive by manual quality check on the Interactive Genome Viewer, and they were therefore removed from downstream analysis. For each of the remaining 10 genes, the number of case carriers ranged from 5 to 10 (3.6%–7.19%) (Table 3). The effect sizes estimated by Fisher’s exact tests were large with low precision due to small sample size (Figure 2). These 10 genes were analysed further to evaluate their associations with AFFs.

Table 3 Discovery and replication results for the 10 suggested genes in hypothesis free analysis.

	Study	# Carriers in AFF cases (%)	# Carriers in controls with BP use (%)	# Carriers in population controls (%)	Compared to controls with BP use	Compared to population controls	
					SKAT-O p-valuesa	ORs (95 % CI)b	SKAT-O p-valuesa	ORs (95 % CI)b	
CYB5D1	Discovery
Replication	5 (3.60%)
0 (0.00%)	0 (0.00%)
0 (0.00%)	20 (0.77%)
6 (0.77%)	0.0006
0.06	Inf (1.31-Inf)
NA	0.005
0.7	3.80 (0.93-11.55)
0.00 (0.00-12.73)	
XRN2	Discovery
Replication	10 (7.19%)
3 (5.66%)	1 (0.51%)
1 (2.56%)	61 (2.36%)
30 (3.84%)	0.002
0.6	15.02 (2.09-657.20)
2.26 (0.17-122.71)	0.0007
0.03	3.21 (1.43-6.49)
1.50 (0.28-5.09)	
SORD	Discovery
Replication	10 (7.19%)
1 (1.89%)	2 (1.02%)
0 (0.00%)	37 (1.43%)
16 (2.05%)	0.004
0.07	7.48 (1.56-71.45)
Inf (0.14-Inf)	0.0000000006
0.8	5.33 (2.31-11.24)
0.92 (0.02-6.15)	
CCDC22	Discovery
Replication	6 (4.32%)
3 (5.66%)	1 (0.51%)
6 (15.38%)	21 (0.81%)
35 (4.48%)	0.008
0.9	8.75 (1.04-405.85)
0.33 (0.05-1.69)	0.01
1	4.55 (1.32-12.66)
1.28 (0.24-4.28)	
DDX4	Discovery
Replication	7 (5.04%)
2 (3.77%)	2 (1.02%)
2 (5.13%)	40 (1.55%)
28 (3.59%)	0.004
0.1	5.12 (0.96-51.25)
0.73 (0.05-10.47)	0.8
1	0.93 (0.11-3.65)
1.05 (0.12-4.39)	
ZNF773	Discovery
Replication	6 (4.32%)
0 (0.00%)	0 (0.00%)
1 (2.56%)	48 (1.85%)
8 (1.02%)	0.0008
0.1	Inf (1.70-Inf)
0.00 (0.00-28.70)	0.05
0.6	2.38 (0.82-5.72)
0.00 (0.00-8.78)	
BRIP1	Discovery
Replication	6 (4.32%)
1 (1.89%)	0 (0.00%)
2 (5.13%)	114 (4.40%)
34 (4.35%)	0.003
0.6	Inf (1.70-Inf)
0.36 (0.01-7.15)	0.6
0.7	0.98 (0.35-2.25)
0.42 (0.01-2.63)	
GBP4	Discovery
Replication	7 (5.04%)
0 (0.00%)	1 (0.51%)
1 (2.56%)	93 (3.59%)
19 (2.43%)	0.003
0.6	10.28 (1.30-467.32)
0.00 (0.00-28.70)	0.2
0.4	1.42 (0.54-3.13)
0.00 (0.00-3.19)	
PHAX	Discovery
Replication	8 (5.76%)
0 (0.00%)	0 (0.00%)
1 (2.56%)	25 (0.97%)
16 (2.05%)	0.004
0.2	Inf (2.49-Inf)
0.00 (0.00-28.70)	0.0001
0.4	6.24 (2.39-14.63)
0.00 (0.00-3.87)	
PLOD2	Discovery
Replication	6 (4.32%)
4 (7.55%)	0 (0.00%)
1 (2.56%)	51 (1.97%)
17 (2.18%)	0.0006
0.06	Inf (1.70-Inf)
3.07 (0.29-156.60)	0.1
0.09	1.85 (0.57-4.72)
3.66 (0.86-11.82)	
a Corrected for sex and four principal components;

b Results by Fisher’s exact test.

Figure 2 Odds ratios and 95% CI for the 10 suggested genes estimated by Fisher’s exact test, results for cases compared with 196 BP-use matched controls and 2588 unmatched controls, respectively.

Downstream analyses for the 10 suggested genes

The 10 genes with suggestive signals in the discovery dataset were investigated using the subsequent downstream methods, and their evidence was combined to determine the likelihood of these genes being associated with AFFs. These genes were prioritized accordingly.

Compared with unmatched controls

Five out of ten genes showed suggestive signals (p-value <.01) when investigated by comparing the same 139 European cases to 2588 unmatched European controls (including the 196 matched controls), namely CYB5D1, XRN2, SORD, CCDC22, and PHAX. The estimated ORs are shown in Figure 2.

Carriers in non-European cases

Among the 46 non-European cases, 22 were of Asian origin. As shown in Figure 3, for SORD, BRIP1, PLOD2, GBP4, DDX4, and ZNF773, a similar or higher percentage of carriers of rare, protein-altering variants were observed in both Asian cases and cases of other ethnicities, compared with European cases. Specifically, for PLOD2, BRIP1, and SORD, a much higher percentage of carriers of rare, protein-altering variants were observed in Asian cases compared with European cases. Lower percentages were observed for XRN2 and PHAX in the Asian cases compared with European cases, but their frequencies were still higher than the frequencies in matched and unmatched European controls (Figure 3). For gene CCDC22 and CYB5D1, no carriers were present in Asian cases or cases of other ethnicities.

Figure 3 Comparison of percentage of carriers in 139 European cases, 22 Asian cases, and 24 cases of other ethnicities, 196 European controls matched for bisphosphonate (BP) use, and 2588 European controls unmatched for BP use.

CADD scores of variants in cases and controls

The number of variants in the cases left after filtering was all below 5 for the 10 suggested genes, as shown in the number of points displayed in Figure 4. Each variant was carried by one or more case. Seven out of ten genes had variants with lower or similar CADD scores in cases compared with those in matched or unmatched controls. In contrast, CYB5D1, XRN2, and CCDC22 showed higher CADD scores for variants in cases than in controls.

Figure 4 CADD scores of variants in cases compared with matched and unmatched controls, respectively.

Comparison with previous genetic findings

None of the 10 suggested genes were among the 132 genes with pathogenic variants present in AFF patients previously reported by Garcia-Giralt et al.32 When these 10 genes were compared with the results of previous genetic studies of AFFs, we found only one gene, SORD, that was linked to a previous finding (Table 4). This gene was within 1 Mb from rs62026663, a suggestive signal identified by the previous GWAS of AFF.18 Three SNPs associated with SORD expression in three different tissues, according to the GTEx Portal, were in linkage disequilibrium (LD; R2 ≥ 0.29) with rs62026663 (known as significant eQTLs) (Table 5).

Gene annotation of suggested genes

Annotations of the 10 suggested genes are shown in Table S3. Seven genes, that is, CCDC22, CYB5D1, GBP4, PHAX, PLOD2, XRN2, and ZNF773, showed expression both in osteoblasts or osteoclasts, and in the Osteogene samples. The expression of PHAX and XRN2 were negatively associated with the expression of haemoglobins (r < −0.6) in Osteogene samples, which contain both bone tissues and blood. Since haemoglobins are mainly expressed in red blood cells, this suggests that PHAX and XRN2 are likely to be specifically expressed in bone tissues. PLOD2 and DDX4 were annotated to GO terms that were over-represented for genes associated with monogenic bone disorders (GO:0007275 multicellular organism development and GO:0030198 extracellular matrix organization, respectively). CYB5D1 was associated with decreased bone mineral content in the IMPC database.

Replication in a Swedish cohort

The 10 genes with suggestive signals in the discovery dataset were analysed using the same methods in Swedish dataset comprising 53 AFF cases and 781 unmatched controls, with a subset analysis of 39 controls matched for a diagnosis of osteoporosis or BP/denosumab use. Four genes had more than one carrier of rare, protein-altering variants among cases, but only XRN2 and PLOD2 showed a higher prevalence of these variants in cases when compared with both matched and unmatched controls. None of the 10 suggested genes in the hypothesis free-approach showed a tendency to an association with AFF in the replication dataset (p-value >.01, shown in Table 3).

Prioritization of suggested genes

Based on a review of all the results given above using scores as presented in Table S4, XRN2, SORD2, and PLOD2 ranked among the most likely novel candidate genes associated with AFF, followed by CYB5D1, CCDC22, and PHAX.

Discussion

In the current study, we analysed genetic data from the largest genetic dataset for AFF so far, to our knowledge: 139 European cases and 46 non-European cases in the discovery dataset and 53 European cases in the replication dataset, totaling 248 AFF cases. We investigated, for the first time, the contributions of rare genetic variants for AFFs using a gene-based approach. In the initial analysis, we compared 139 European cases with 196 controls matched for BP use for validation of known candidate genes and exploration of novel genes. Finally, based on additional analyses, we prioritized three of these genes as the most likely candidates associated with AFF, namely XRN2, SORD, and PLOD2.

Candidate gene analysis

In an analysis of 203 candidate genes, only PLOD2 showed suggestively significant results when analysing rare, protein-altering variants with CADD scores greater than 0 and 10, but not when analysing variants with CADD greater than 20. In our previous study, several genetic variants associated with monogenic bone disorders were identified, such as in CTSK, PLS3, COL1A1, COL1A2, LRP5, and ALPL.7,11 In this study, carriers of these genetic variants were not sufficiently enriched in the AFF cases compared with controls. With a limited number of carriers for each of these genes (Table S2), the current study lacks sufficient power to evaluate these genes. No evidence of associations with AFF was observed for other genes associated with monogenic bone disorders, nor those associated with bone mineralization, such as ENPP1, SLC34A1, and SLC9A3R1, as previously suggested by Furukawa et al. and Marini et al.12,13 We found no association between AFF and genes involved in the mechanism of BP action, such as GGPS1 and CYP1A1, which were suggested in the study of three AFF sisters by Roca-Ayats et al., nor with ATRAID, in which a variant was indicated to be related to AFF and increased sensitivity to BPs.14,16 While multiple genes have been implicated in AFFs, our results show that previous genetic findings of AFFs are not widely shared among most AFF cases, underscoring that genetic background of AFF is heterogenous.

Table 4 SORD was lined to a previous genetic study of AFFs.

Gene	Current study	Previous genetic finding	
	# Case carriers (%)	# Control carriers (%)	Study	Finding	Remark	
SORD	10 (6.9%)	2 (1%)	GWAS of AFF by Kharazmi and colleagues18	A suggestive SNP rs62026663 is in LD with three other SNPs (R2 ≥ 0.3) that are associated with SORD expression in three different tissues (GTEx)	14 suggestive SNPs (p-value <1e−5) were identified in the GWAS of AFF	

Table 5 SNPs in LD with rs62026663 and in significant eQTL with the SORD gene

SNP in LD	R 2  with rs62026663	p-value	NES	Tissue	
rs56351200	0.29	2.20E−07	0.12	Cells - Cultured fibroblasts	
rs62025132	0.33	6.10E−06	0.22	Artery - Aorta	
rs62025168	0.29	2.30E−05	0.20	Skin - Not Sun Exposed (Suprapubic)	
SNP, single nucleotide polymorphism; eQTL, expression quantitative trait loci; LD, linkage disequilibrium; NES, normalized effect size.

Hypothesis-free analysis approach

In the hypothesis-free approach, we further evaluated the 10 suggested genes although no genome-wide significant results were identified after correction for multiple testing. The results were subject to an inflated type I error rate due to multiple testing which may yield false positive findings. However, if rare, protein-altering variants in a gene are associated with AFF, it is likely that the gene would have emerged as a suggested gene in our analysis. It should be noted that the current study lacked sufficient power to produce statistically significant outcomes.

Therefore, several methods were used to assist finding the true associations among the suggested genes. These genes were (1) investigated by comparing cases with unmatched controls; (2) analysed in an independent replication dataset, (3) compared with previous genetic findings, (4) assessed in non-European samples, including a subgroup of Asian cases, and (5) evaluated with CADD scores and functional annotation. Finally, we prioritized the genes by integrating the methods. Genes XRN2, SORD, and PLOD2 were considered the most likely candidates associated with AFF and are discussed below for their potential involvement in AFF pathogenesis. The genes CYB5D1, CCDC22, and PHAX have not been reported in the literature to have any known involvement in bone.

The gene XRN2 encodes the enzyme 5'-3' exoribonuclease 2, which may promote transcription termination by RNA polymerase II.40 We found 10 carriers (7%) among cases in the discovery dataset. In contrast, only 2% carrier rate was observed in the unmatched controls in the discovery dataset. The variants carried by cases tend to be more damaging (higher CADD scores) than variants carried by controls. In the replication dataset, we also observed higher prevalence of carriers of XRN2 variants in cases than in controls (6% in cases and 4% in unmatched controls). The XRN2 gene is expressed in both human osteoblasts and osteoclasts, and specifically expressed in bone tissues in the Osteogene database, suggesting an important function in bone metabolism. It has been suggested that XRN2 can bind to NF-κB repressing factor (NKRF),41 a protein which is described as PLS3 interactor. NKRF was found to increasingly translocate to the nucleus upon overexpression of PLS3.42,43 Previously, AFF has been reported in patients with X-linked osteoporosis caused by PLS3 variants.11 Although the role of PLS3 variants in AFFs has not be established, it is possible that a pathogenic XRN2 variant may result in AFF through the same pathway as the effect of PLS3 variants.

The gene SORD encodes sorbitol dehydrogenase, an enzyme that converts sorbitol to fructose.44 We found SORD variants in 10 AFF cases (7%) in the discovery dataset but only one AFF case (2%) in the replication dataset. We also observed a higher proportion of case carriers of SORD variants in the Asian cases (13%) than in European cases. The gene was linked to an SNP (rs62026663) suggestively associated with AFF identified by a previous GWAS through eQTL and SNPs in LD.18 Patients with biallelic SORD variants had sorbitol dehydrogenase deficiency with peripheral neuropathy that mainly affects lower limbs.45 The condition has been associated with accumulated intracellular and serum sorbitol levels.45 Increased sorbitol has been suggested to elevate advanced glycation end products (AGEs) in serum and tissues such as lenses and liver.44,46 AGE cross-linking in bone has been associated with higher fracture risk,47 suggesting that SORD may play a role in the maintenance of bone quality. These effects could be compounded by prolonged BP use, which has also been associated with increased AGE accumulation in animal bone.48

Biallelic pathogenic variants in the gene PLOD2 are responsible for the autosomal recessive Bruck syndrome 2, which is a form of osteogenesis imperfecta caused by a deficiency of telopeptide lysyl hydroxylase in bone collagen.49 This gene was included in the candidate gene analysis for its association with a monogenic bone disorder. We observed 4.3% and 7.5% of cases carrying variants in PLOD2 in the discovery dataset and the replication dataset, respectively. In contrast, only 2% in the unmatched controls carried PLOD2 variants in both datasets. Notably, no significant trend of ORs was observed when the gene was analysed in steps based on CADD score thresholds at 0, 10, and 20, and the gene only showed suggestive signals in the first two steps. However, PLOD2 remains a plausible candidate for AFF. PLOD2 rare variants were particularly more prevalent in the Asian AFF cases. One of the Asian cases, previously reported by our group, was a South-East Asian AFF patient whose mother also had AFF but was not genetically tested.10

Proposal of genetic architecture of AFFs

AFF cases may have a monogenic background, as AFF has been associated with monogenic bone disorders, and in some cases appears to cluster in families.7,11 There is also a high prevalence of suspected monogenic bone disorders in AFF patients.11 Moreover, it is possible that additional genes may harbor monogenic causes of AFF, which have not yet been identified in AFF patients, since the genetic cause for 50% of the AFF patients with suspected monogenic bone disorder remains unknown/uncertain.11 However, our results showed that there was no sufficient enrichment of potentially pathogenic variants in any of the genes analysed in our cases compared with controls. These findings suggest that sporadic AFF cases have a high level of genetic heterogeneity, more than that observed in monogenic bone disorders, as a single genetic factor responsible for AFFs is not commonly shared among AFF cases. Moreover, in some cases, AFF could also have an oligogenic background, being caused by the combined or interacted effects of multiple genes, each with milder effects, and environmental factors. Above all, multiple distinct genotypes, potentially involving genotype-environment interactions, could lead to AFF. Even though we analysed the largest genetic dataset of AFF patients so far, the current study is still relatively small and may not have enough power to detect these associations because of the genetic heterogeneity of the AFF phenotype. Larger genetic studies are needed to understand the genetic architecture of AFF.

AFF patients may have various primary conditions caused by different genetic mechanisms, such as an underlying genetic bone disorder, an underlying disorder for which glucocorticoids are prescribed, an altered response to BP treatment (eg BP sensitivity), a genetic predisposition to a specific type of geometry of the femur, or a combination of these factors, which result in bone properties allowing susceptibility to AFF. For instance, several studies have suggested that bone biopsies from AFF patients show higher bone mineralization and higher enzymatic collagen maturity as compared with non-AFF patients.50 Further studies are needed to confirm the underlying pathogenic mechanisms for AFF.

Strengths and limitations

In this study, we present the largest genetic study of AFF patients to date and have analysed genetic data by a gene-based association approach. This approach aggregated various possibly relevant rare variants within a gene in different cases and could identify genes with causal variants not detected in single-variant analyses due to limited power. However, this study has limitations. First, despite being the largest to date, the study is still underpowered and is subject to an inflated type I error rate due to low total sample size and genetic heterogeneity. Nevertheless, by integrating multiple downstream analyses for gene prioritization, we showed that a few genes might be true associations with AFFs and worthy of further investigation. It is important to emphasize that these results should be interpreted with caution, and it is necessary to replicate the suggested genes before proceeding to validate them through functional experiments. Secondly, we analysed all the variants within each gene equally without applying weight, although it is expected that rarer and more damaging variants should have a greater effect than less rare, damaging variants. However, it would be challenging to predict the pathogenicity of each variant and estimate the differences between them, especially considering that it could vary across different genes. Thirdly, due to the use of different exome capture kits for cases and controls and the subsequent need to harmonize the datasets, not all possible exonic information could be extracted. Consequently, some relevant variants may have been missed. Moreover, it was assumed that controls were not affected by AFF, which is likely to be true given the rarity of AFF in the population. Furthermore, the study only focused on rare, protein altering variants that are more likely to influence gene function. We may have missed other genetic variants that contribute to AFF, such as common genetic variants and variants in regulatory regions. In future studies, the combined effect of rare and common genetic variants on AFF could be considered. We acknowledge potential biases, including incomplete adjustment for indication of BP use, given that matched controls have shorter BP use duration. Additionally, selection bias may arise as our cohort may not fully represent the broader AFF population, particularly considering that cases collected at bone expertise centers in the Netherlands may represent more severe or complicated AFF cases. Finally, due to limited power, we were unable to further explore interactions between genetic variants and BP use or perform subgroup analyses based on potential contributing factors such as glucocorticoid use.

Conclusion

In summary, gene-based association analyses of 203 candidate genes in 139 AFF cases and 196 BP-matched controls showed no significant enrichment for carriers of rare, protein-altering variants in AFF cases compared with controls. In the-hypothesis free approach, we tested a total of 13 198 genes and identified 10 suggested genes. Although none of the genes showed suggestive signals in a replication dataset of 53 AFF cases, we identified the genes PLOD2 and XRN2 with consistent direction of effect estimates, and the gene SORD by comparing the gene list with previous genetic findings in AFFs. Moreover, PLOD2 and SORD also presented with higher carrier frequencies in 22 Asian cases, again suggestive of a true association. We thus conclude that XRN2, PLOD2, and SORD are potential novel candidates associated with AFF. The results further suggest that AFF may be caused by different genes in different individuals, thus strengthening the notion that genetic heterogeneity may be important. A much larger sample size is now needed to fully understand the genetic architecture of AFFs.

Supplementary Material

Supplementary_figures_Sept_2023_zjae122

Table_S1-S4_APR2024_zjae122

Supplementary_methods_APR2024_zjae122

Acknowledgments

We are grateful to Sandra Smits and Els van Zaanen (Erasmus MC, The Netherlands) for their practical assistance in the study. We are indebted to Fernando Rivadeneira and Jeroen G.J. van de Peppel (Erasmus MC, the Netherlands) for providing gene expression datasets for our analysis. We thank Joost Verlouw (Erasmus MC, the Netherlands) for his bioinformatic input in the study. We thank Linda Broer (Erasmus MC, the Netherlands) for sharing her scripts for ethnicity annotation. We thank Frances Milat (Monash Health, Australia), Christopher Yates (Western Health, Australia), and Cherie Chiang (Austin Health, Australia) for referring participants to the Australian dataset. Our sincere gratitude extends to all the participants. The generation and management of the exome sequencing data for the RS were executed by the Genomics Core Facility, Erasmus MC, Rotterdam, The Netherlands. We thank Pascal Arp, Mila Jhamai, Jeroen G.J. van Rooij, Marijn Verkerk, and Robert Kraaij for their help in creating the RS-Exome Sequencing database.

It was performed at the single nucleotide polymorphism and sequencing (SNP&SEQ) Technology Platform at Uppsala University that is part of the National Genomics Infrastructure supported by the Swedish Research Council and the Knut and Alice Wallenberg Foundation. Computations were done on resources provided by the Swedish National Infrastructure for Computing through the Uppsala Multidisciplinary Centre for Advanced Computational Science (UPPMAX). The authors are grateful to the study participants, the Swedegene staff, and the senior consultants in orthopaedics that adjudicated the cases.

Author contributions

Wei Zhou (Conceptualization, Data curation, Formal analysis, Methodology, Visualization, Writing—original draft, Writing—review & editing), Joel Ås (Data curation, Formal analysis, Writing—review & editing), Catherine Shore-Lorenti (Data curation, Project administration, Resources, Writing—review & editing), Hanh H. Nguyen (Data curation, Project administration, Resources, Writing—review & editing), Denise M. van de Laarschot (Conceptualization, Project administration, Writing—review & editing), Shoshana Sztal-Mazer (Resources, Writing—review & editing), Vivian Grill (Resources, Writing—review & editing), Christian M. Girgis (Resources, Writing—review & editing), Bruno H. Ch. Stricker (Data curation, Resources, Writing—review & editing), Bram C.J. van der Eerden (Writing—review & editing), Rajesh V. Thakker (Writing—review & editing), Natasha M. Appelman-Dijkstra (Resources, Writing—review & editing), Mia Wadelius (Resources, Supervision, Writing—review & editing), Roderick J. Clifton-Bligh (Resources, Writing—review & editing), Pär Hallberg (Resources, Supervision, Writing—review & editing), Annemieke J.M.H. Verkerk (Project administration, Supervision, Writing—review & editing), Jeroen G.J. van Rooij (Conceptualization, Methodology, Supervision, Writing—review & editing), Peter R. Ebeling (Conceptualization, Funding acquisition, Investigation, Resources, Writing—review & editing), and M. Carola Zillikens (Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Writing—review & editing)

Funding

Jaap Schouten Foundation, Rotterdam, The Netherlands, and National Health and Medical Research Council of Australia (GNT1197958) provided funding for this study. R.V.T. has received a Wellcome Trust Investigator Award (grant number 106995/Z/15/Z), National Institute for Health Research (NIHR) Oxford Biomedical Research Centre Programme, and NIHR Senior Investigator Award (grant number NF- SI-0514–10091). Generation of the Exome Sequencing data set of the Rotterdam Study was funded by the Netherlands Genomics Initiative (NGI)/Netherlands Organisation for Scientific Research (NWO) sponsored Netherlands Consortium for Healthy Aging (NCHA; project nr. 050-060-810), by the Genetic Laboratory of the Department of Internal Medicine, Erasmus MC, and by a Complementation Project of the Biobanking and Biomolecular Research Infrastructure Netherlands (BBMRI-NL; www.bbmri.nl; project number CP2010-41). The Rotterdam Study is funded by Erasmus Medical Center and Erasmus University, Rotterdam, Netherlands Organization for the Health Research and Development (ZonMw), the Research Institute for Diseases in the Elderly (RIDE), the Ministry of Education, Culture and Science, the Ministry for Health, Welfare and Sports, the European Commission (DG XII), and the Municipality of Rotterdam. Collection of Swedegene’s replication cohort was supported by the Swedish Research Council (Medicine 521-2011-2440, 521-2014-3370, and 2018-03307) and Clinical Research Support (Avtal om Läkarutbildning och Forskning, ALF) at Uppsala University, Sweden. Whole genome sequencing was granted by the Science for Life Laboratory’s Swedish Genomes Program 2017 that was supported by the Knut and Alice Wallenberg Foundation (application ID NP:00085). The funding agencies played no role in the study design or in data collection and analyses.

Conflicts of interest

PRE: Research funding from Amgen, Sanofi, and Alexion. Honoraria from Amgen, Pfizer and Gedeon Richter. NMA-D: Lecture fees from Amgen, UCB. Research funding Takeda, Kyowa Kirin. RVT: has received grants from Novartis Pharma AG, Novo Nordisk, and the Marshall Smith Syndrome Foundation for unrelated studies. RCB: lecture fees and/or advisory board fees from Amgen, Eisai, Ipsen, Kyowa Kirin.

Data availability

Sharing raw or processed individualized sequencing results of the patients are not allowed due to General Data Protection Regulation (GDPR). Requests to access the datasets should be directed to MCZ, m.c.zillikens@erasmusmc.nl.
==== Refs
References

1. Jha  S, Wang  Z, Laucis  N, Bhattacharyya  T. Trends in media reports, oral bisphosphonate prescriptions, and hip fractures 1996–2012: an ecological analysis. J Bone Miner Res. 2015;30 (12 ):2179–2187. 10.1002/jbmr.2565 26018247
2. Wysowski  DK, Greene  P. Trends in osteoporosis treatment with oral and intravenous bisphosphonates in the United States, 2002–2012. Bone. 2013;57 (2 ):423–428. 10.1016/j.bone.2013.09.008 24063946
3. Khow  KSF, Shibu  P, Yu  SCY, Chehade  MJ, Visvanathan  R. Epidemiology and postoperative outcomes of atypical femoral fractures in older adults: a systematic review. J Nutr Health Aging. 2017;21 (1 ):83–91. 10.1007/s12603-015-0652-3 27999854
4. Black  DM, Geiger  EJ, Eastell  R, et al.  Atypical femur fracture risk versus fragility fracture prevention with bisphosphonates. N Engl J Med. 2020;383 (8 ):743–753. 10.1056/NEJMoa1916525 32813950
5. Shane  E, Burr  D, Abrahamsen  B, et al.  Atypical subtrochanteric and diaphyseal femoral fractures: second report of a task force of the American Society for Bone and Mineral Research. J Bone Miner Res. 2014;29 (1 ):1–23. 10.1002/jbmr.1998 23712442
6. Soo-Young  J, Kyu-Tae  H, Chang-Wug  O, et al.  Mid-term outcomes after the surgical treatment of atypical femoral fractures. Bone Joint J. 2021;103-B (11 ):1648–1655. 10.1302/0301-620X.103B11.BJJ-2021-0416.R1 34719278
7. Zhou  W, van Rooij  JGJ, Ebeling  PR, Verkerk  AJMH, Zillikens  MC. The genetics of atypical femur fractures—a systematic review. Curr Osteoporos Rep. 2021;19 (2 ):123–130. 10.1007/s11914-021-00658-y 33587247
8. Nguyen  HH, van de Laarschot  DM, Verkerk  AJ, Milat  F, Zillikens  MC, Ebeling  PR. Genetic risk factors for atypical femoral fractures (AFFs): a systematic review. JBMR Plus. 2018;2 (1 ):1–11. 10.1002/jbm4.10024 30283886
9. Roca-Ayats  N, Balcells  S, Garcia-Giralt  N, et al.  GGPS1 mutation and atypical femoral fractures with bisphosphonates. N Engl J Med. 2017;376 (18 ):1794–1795. 10.1056/NEJMc1612804 28467865
10. Zhou  W, Nguyen  HH, van de Laarschot  DM, et al.  Whole exome sequencing in two southeast Asian families with atypical femur fractures. JBMR Plus. 2022;6 (8 ):e10659. 10.1002/jbm4.10659 35991532
11. Zhou  W, van Rooij  JG, van de Laarschot  DM, et al.  Prevalence of monogenic bone disorders in a Dutch cohort of atypical femur fracture patients. J Bone Miner Res. 2023;38 (6):896–906. 10.1002/jbmr.4801
12. Furukawa  H, Oka  S, Kondo  N, et al.  The contribution of deleterious rare alleles in ENPP1 and osteomalacia causative genes to atypical femoral fracture. J Clin Endocrinol Metab. 2022;107 (5 ):e1890–e1898. 10.1210/clinem/dgac022 35038731
13. Marini  F, Giusti  F, Marasco  E, et al.  High frequency of heterozygous rare variants of the SLC34A1 and SLC9A3R1 genes in patients with atypical femur fracture. Eur J Endocrinol. 2023;188 (1 ):176–188. 10.1093/ejendo/lvad001
14. Roca-Ayats  N, Ng  PY, Garcia-Giralt  N, et al.  Functional characterization of a GGPPS variant identified in atypical femoral fracture patients and delineation of the role of GGPPS in bone-relevant cell types. J Bone Miner Res. 2018;33 (12 ):2091–2098. 10.1002/jbmr.3580 30184270
15. Ugartondo  N, Martínez-Gil  N, Esteve  M, et al.  Functional analyses of four CYP1A1 missense mutations present in patients with atypical femoral fractures. Int J Mol Sci. 2021;22 (14 ):7395. 10.3390/ijms22147395
16. Surface  LE, Burrow  DT, Li  J, et al.  ATRAID regulates the action of nitrogen-containing bisphosphonates on bone. Sci Transl Med. 2020;12 (544 ):eaav9166. 10.1126/scitranslmed.aav9166
17. Pérez-Núñez  I, Pérez-Castrillón  JL, Zarrabeitia  MT, et al.  Exon array analysis reveals genetic heterogeneity in atypical femoral fractures. A pilot study Mol Cell Biochem. 2015;409 (1-2 ):45–50. 10.1007/s11010-015-2510-3 26160281
18. Kharazmi  M, Michaëlsson  K, Schilcher  J, et al.  A genome-wide association study of bisphosphonate-associated atypical femoral fracture. Calcif Tissue Int. 2019;105 (1 ):51–67. 10.1007/s00223-019-00546-9 31006051
19. Li  B, Leal  SM. Methods for detecting associations with rare variants for common diseases: application to analysis of sequence data. Am J Hum Genet. 2008;83 (3 ):311–321. 10.1016/j.ajhg.2008.06.024 18691683
20. Ikram  MA, Brusselle  G, Ghanbari  M, et al.  Objectives, design and main findings until 2020 from the Rotterdam study. Eur J Epidemiol. 2020;35 (5 ):483–517. 10.1007/s10654-020-00640-5 32367290
21. van Rooij  JGJ, Jhamai  M, Arp  PP, et al.  Population-specific genetic variation in large sequencing data sets: why more data is still better. Eur J Hum Genet. 2017;25 (10 ):1173–1175. 10.1038/ejhg.2017.110 28905877
22. Raghavan  NS, Brickman  AM, Andrews  H, et al.  Whole-exome sequencing in 20,197 persons for rare variants in Alzheimer's disease. Ann Clin Transl Neurol. 2018;5 (7 ):832–842. 10.1002/acn3.582 30009200
23. Wang  Q, Pierce-Hoffman  E, Cummings  BB, et al.  Landscape of multi-nucleotide variants in 125,748 human exomes and 15,708 genomes. Nat Commun. 2020;11 (1 ):2539. 10.1038/s41467-019-12438-5 32461613
24. Rentzsch  P, Witten  D, Cooper  GM, Shendure  J, Kircher  M. CADD: predicting the deleteriousness of variants throughout the human genome. Nucleic Acids Res. 2018;47 (D1 ):D886–D894. 10.1093/nar/gky1016
25. Maróti  Z, Boldogkői  Z, Tombácz  D, Snyder  M, Kalmár  T. Evaluation of whole exome sequencing as an alternative to BeadChip and whole genome sequencing in human population genetic analysis. BMC Genomics. 2018;19 (1 ):778. 10.1186/s12864-018-5168-x 30373510
26. Zhou  W, Nielsen  JB, Fritsche  LG, et al.  Efficiently controlling for case-control imbalance and sample relatedness in large-scale genetic association studies. Nat Genet. 2018;50 (9 ):1335–1341. 10.1038/s41588-018-0184-y 30104761
27. Lee  S, Fuchsberger  C, Kim  S, Scott  L. An efficient resampling method for calibrating single and gene-based rare variant association analysis in case–control studies. Biostatistics. 2015;17 (1 ):1–15. 10.1093/biostatistics/kxv033 26363037
28. Bouaziz  M, Mullaert  J, Bigio  B, et al.  Controlling for human population stratification in rare variant association studies. Sci Rep. 2021;11 (1 ):19015. 10.1038/s41598-021-98370-5 34561511
29. Robinson  JT, Thorvaldsdóttir  H, Wenger  AM, Zehir  A, Mesirov  JP. Variant review with the integrative genomics viewer. Cancer Res. 2017;77 (21 ):e31–e34. 10.1158/0008-5472.CAN-17-0337 29092934
30. Mol  MO, van der Lee  SJ, Hulsman  M, et al.  Mapping the genetic landscape of early-onset Alzheimer’s disease in a cohort of 36 families. Alzheimers Res Ther. 2022;14 (1 ):77. 10.1186/s13195-022-01018-3 35650585
31. Kircher  M, Witten  DM, Jain  P, O'Roak  BJ, Cooper  GM, Shendure  J. A general framework for estimating the relative pathogenicity of human genetic variants. Nat Genet. 2014;46 (3 ):310–315. 10.1038/ng.2892 24487276
32. Garcia-Giralt  N, Roca-Ayats  N, Abril  JF, et al.  Gene network of susceptibility to atypical femoral fractures related to bisphosphonate treatment. Genes (Basel). 2022;13 (1 ):146. 10.3390/genes13010146
33. Prijatelj  V, Reppe  S, Dietz  M, et al.  Differential expression analysis of osteoarthritic femoral head bone fragments uncovers underlying transcriptomic determinants. Bone Reports. 2020;13 :100503. 10.1016/j.bonr.2020.100503
34. Hallberg  P, Yue  QY, Eliasson  E, Melhus  H, Ås  J, Wadelius  M. SWEDEGENE-a Swedish nation-wide DNA sample collection for pharmacogenomic studies of serious adverse drug reactions. Pharmacogenomics J. 2020;20 (4 ):579–585. 10.1038/s41397-020-0148-3 31949290
35. Li  H, Durbin  R. Fast and accurate short read alignment with burrows-wheeler transform. Bioinformatics. 2009;25 (14 ):1754–1760. 10.1093/bioinformatics/btp324 19451168
36. Li  H, Handsaker  B, Wysoker  A, et al.  The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25 (16 ):2078–2079. 10.1093/bioinformatics/btp352 19505943
37. McKenna  A, Hanna  M, Banks  E, et al.  The genome analysis toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20 (9 ):1297–1303. 10.1101/gr.107524.110 20644199
38. DePristo  MA, Banks  E, Poplin  R, et al.  A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat Genet. 2011;43 (5 ):491–498. 10.1038/ng.806 21478889
39. Poplin  R, Ruano-Rubio  V, DePristo  MA, et al.  Scaling accurate genetic variant discovery to tens of thousands of samples. bioRxiv. 2018;201178. 10.1101/201178
40. Skourti-Stathaki  K, Proudfoot  NJ, Gromak  N. Human senataxin resolves RNA/DNA hybrids formed at transcriptional pause sites to promote Xrn2-dependent termination. Mol Cell. 2011;42 (6 ):794–805. 10.1016/j.molcel.2011.04.026 21700224
41. Alexandrova  J, Piñeiro  D, Jukes-Jones  R, Mordue  R, Stoneley  M, Willis  AE. Full-length NF-κB repressing factor contains an XRN2 binding domain. Biochem J. 2020;477 (4 ):773–786. 10.1042/BCJ20190733 32011671
42. Besio  R, Chow  C-W, Tonelli  F, Marini  JC, Forlino  A. Bone biology: insights from osteogenesis imperfecta and related rare fragility syndromes. FEBS J. 2019;286 (15 ):3033–3056. 10.1111/febs.14963 31220415
43. Neugebauer  J, Heilig  J, Hosseinibarkooie  S, et al.  Plastin 3 influences bone homeostasis through regulation of osteoclast activity. Hum Mol Genet. 2018;27 (24 ):4249–4262. 10.1093/hmg/ddy318 30204862
44. Moemen  LA, Abdel Hamid  MA, Wahab  SA, et al.  Role of advanced glycation end products and sorbitol dehydrogenase in the pathogenesis of diabetic retinopathy. Bull Natl Res Cent. 2020;44 (1 ):58. 10.1186/s42269-020-00304-0
45. Cortese  A, Zhu  Y, Rebelo  AP, et al.  Biallelic mutations in SORD cause a common and potentially treatable hereditary neuropathy with implications for diabetes. Nat Genet. 2020;52 (5 ):473–481. 10.1038/s41588-020-0615-4 32367058
46. Gugliucci  A . Formation of fructose-mediated advanced glycation end products and their roles in metabolic and inflammatory diseases. Adv Nutr. 2017;8 (1 ):54–62. 10.3945/an.116.013912 28096127
47. Suzuki  A, Yabu  A, Nakamura  H. Advanced glycation end products in musculoskeletal system and disorders. Methods. 2022;203 :179–186. 10.1016/j.ymeth.2020.09.012 32987130
48. Saito  M, Mori  S, Mashiba  T, Komatsubara  S, Marumo  K. Collagen maturity, glycation induced-pentosidine, and mineralization are increased following 3-year treatment with incadronate in dogs. Osteoporos Int. 2008;19 (9 ):1343–1354. 10.1007/s00198-008-0585-3 18373056
49. van der Slot  AJ, Zuurmond  AM, Bardoel  AF, et al.  Identification of PLOD2 as telopeptide lysyl hydroxylase, an important enzyme in fibrosis. J Biol Chem. 2003;278 (42 ):40967–40972. 10.1074/jbc.M307380200 12881513
50. Ural  A . Biomechanical mechanisms of atypical femoral fracture. J Mech Behav Biomed Mater. 2021;124 :104803. 10.1016/j.jmbbm.2021.104803 34479108
