This article is mentioned in:

Introduction

The incidence of fractures is highest among the young and older populations (1). Of note, fractures in pediatric patients are common injuries that result in pain, as well as short- and long-term morbidity; it has been estimated that 27-50% of all individuals experience a fracture prior to turning 18 years old (2,3). Goulding et al (4) estimated the risk of recurrent fracture to be 12-20% in pediatrics patients aged <16 years. Furthermore, these figures may indicate that pediatric patients who do not have any diseases affecting bone integrity may still be at risk to suffer a bone fracture or even recurrent fractures (2-4). The causes for this may arise from the diet, preferred free-time activities or sex, but may also be linked to genetic factors (3).

The causes and etiopathogenetic pathways of pediatric fractures are not well studied (5). It has been previously reported that bone-associated and bone-independent factors contribute to the fracture risk (2,6). Furthermore, bone mineral content (BMC) and bone mineral density (BMD) have been indicated to be highly heritable in genetic studies (7-11). Mora and Gilsanz (12) revealed that 80% of the variance of the BMC is influenced by genetics. Furthermore, BMD has been reported to be a risk factor for fractures in adults (1), but there is limited research on pediatric fractures. Previous studies have linked a low BMD to increased fracture risk during childhood (3,13); however, in other studies, this association was not identified (6,14). In adults, genome-wide association studies (GWAS) have identified genes from ≥11 different chromosomes (chrs) that may affect BMD either locally or throughout the body (15-21). However, with regards to pediatrics, only few studies have reported on genetic loci that influence BMD or cortical bone thickness (7,8). Furthermore, the heritability of the fractures and BMD differ, and thus, it has been suggested that fracture risk should be analyzed independently of BMD (22).

BMD and BMC are useful when examining the factors underlying fractures (23,24); however, their association with fractures in childhood remains to be fully elucidated. Therefore, the present study focused on fractures rather than BMD or BMC. To the best of our knowledge, only few studies have examined the association between genetics and childhood fracture risk. Therefore, the aim of the present study was to identify whether any genetic loci are associated with fractures prior to the age of 7 years using GWAS.

Patients and methods

Patient population

The study population was a part of the Northern Finland Birth Cohort 1986 (NFBC1986) that originally included all pregnant females (n=9,362) living in the geographic area of the two northernmost provinces of Finland with expected date of delivery between 1st July 1985 and 30th June 1986. These mothers and their live-born infants (n=9,432) were followed up regularly since pregnancy. Immigration during this period to northern Finland was minimal, and therefore, the cohort was predominantly Finnish (25). The detailed description of the cohort is provided in the study published by Järvelin et al (26) and on the NFBC website (27).

Of the cohort population, 6,721 consented to the use of their information for research purposes. The information for the individuals who did not provide consent was not available. Study subjects with diseases of the bone, including osteogenesis imperfecta, or possible bone-affecting malignancies, were excluded from the original study population (n=3). Informed consent from the parents was provided orally during antenatal clinic visits. At the 16-year follow-up (in 2001), the cohort participants and their parents provided written informed consent for the use of all of their collected data for scientific purposes. The participation was voluntary.

Genotype data

Genotypes were determined for the cohort members who had provided consent and their blood samples were used for this purpose. The genotype data were available for 3,230 individuals, including 48 cases who had a fracture and 3,182 controls without fractures. Genotyping was performed using the Illumina HumanOmniExpressExome v1.2 bead chip (Illumina, Inc.) and imputation was performed using the 1000G phase 3 reference panel (28). After quality control (QC), the study included 3,515,000 single-nucleotide polymorphisms (SNPs) (28). The QC process is described in detail in a previous study (28). The data was controlled based on genotype call rate, minor allele frequency, multidimensional scaling outliers, relatedness, heterozygosity rate and Hardy-Weinberg Equilibrium. The imputation quality threshold was 0.4. The flow chart for inclusion is presented in Fig. 1.

Figure 1 Flow chart of the inclusion criterion.

Fractures

The National Hospital Discharge Register (NHDR) provided the data on the fractures (the name of the register has changed and it is now known as Care Register for Health Care). The register contains information on all in-hospital treated fractures with diagnostic codes and discharge dates. The NHDR is coordinated by the Finnish institute for health and welfare and it contains all hospital admissions in Finland (29). This register is the oldest nationwide discharge register worldwide and thoroughly covers acute injuries of the population (annual coverage of injuries, ~100%) and classifies these injuries accurately (annual accuracy, >95%) (30).

In the study area, pediatric fractures are exclusively treated in public hospitals and it is mandatory to send information of these injuries to the NHDR. During the study period, out-patient visits were not recorded in the NHDR and any patients treated outside the hospital were not included. Simple fractures, including non-displaced fracture of wrist or foot, were not included in the NHDR, as these are treated by general practitioners in primary healthcare centers. All the in hospital treated fractures were recognized for the analysis and the patients were classified according to whether they had a fracture or not. The International Classification of Diseases, version-9 (ICD-9) was used to identify the fractures (31). If the patient had several hospitalizations with the same ICD code, the fractures were determined to be caused by ≥2 separate injuries if there was an interval of ≥6 months between the hospitalization dates. The ICD-9 codes reported for the fractures are presented in Table I.

Table I Reported ICD-9 codes.

Table I

Reported ICD-9 codes.

ICD-9 code	Explanation
800	Fracture of vault of skull
802	Fracture of facial bones
810	Fracture of clavicle
812	Fracture of humerus
813	Fracture of radius or ulna
816	Fracture of one or more phalanges of the hand
821	Fracture of other or unspecified parts of femur
823	Fracture of tibia or fibula

[i] ICD, International Classification of Diseases.

The cause of the fracture may also be included in the NHDR data. Out of the 48 individuals who had a fracture in the population of the present study, 37 had a reported mechanism. Furthermore, 15 patients had low-energy accidents, including falling onto a plane, tripping, slipping or falling from a height of <1 m. For the other fractures, it was not possible to reliably determine the energy involved in the injury.

Data analyses

The association between specific genomic regions and fractures was assessed in the GWAS using the frequentist association test of the SNPTEST analysis software (version 2.5.2; University of Oxford). In the GWAS analysis, a linear regression model was fitted to test for additive effects of SNPs (genotype dosage) adjusting for sex and population stratification using four principal components. Sex and the first four principal components were used as covariates.

The post-GWAS QC included the following steps: Minor allele frequency <0.05, deviation from Hardy-Weinberg equilibrium P-value <1x10-7. The generally accepted limit of P<5x10-8 for genome-wide significance and P<5x10-5 for suggestive evidence of association were used. Under the null hypothesis, the threshold where one false positive is expected per genome scan is called the suggestive level (32). Manhattan and Quantile-Quantile (Q-Q) plots of the results were generated using R (version 3.5.0) and the qqman package (version 0.1.7) (33). The closest coding and non-coding genes in the identified loci were determined using the University of California Santa Cruz Genome Browser (December 2013), which is a graphical viewing tool of the human genome assembly. The genes were visualized using the LocusZoom visualization software (version 1.4; University of Michigan).

After GWAS, the lead SNPs were analyzed using the public databases GWAS database (GWASdb) (34), GWAS catalog (35), The Genotype-Tissue Expression (GTEx) Portal (36) and RegulomeDB (37). GWASdb and GWAS catalog are databases of genetic variants identified in published GWASs. GTEx Portal is database of tissue-specific gene expression. RegulomeDB is used to annotate variants to known and predicted regulatory elements. Database searches were performed for the locus surrounding the lead SNP (chr10:11491165-12291165). Associated loci were compared with summary statistics from two studies available on the Genetic Factors for Osteoporosis Consortium (GEFOS) website (38,39). Trajanoska et al (38) included 25 cohorts and had discovery dataset of 37,857 fracture cases and 227,116 controls composed. All the subjects were adults and predominantly of European descent (38). Medina-Gomez et al (39) was meta-analysis of total body bone mineral density and included 66,628 individuals, both children and adults of mostly European descent (39).

Results

Fracture demographics

In total, 48 patients had a fracture prior to the age of 7 and one of these patients had two fractures. The characteristics of the study population are presented in Table II. Of the 48 subjects with fractures, 16 (33.3%) were female and 32 (66.7%) were male. The average age at the time of fracture was 4.05 years (range 0-6 years). Furthermore, the fractures included one cervical spine fracture (ICD-9 806), two skull fractures (ICD-9 800), eight lower limb fractures (ICD-9 821, 823) and 37 upper limb fractures (ICD-9 810, 812, 813 and 816).

Table II Characteristics of the study population.

Table II

Characteristics of the study population.

Characteristic	Value
Sex (female/male)
Fracture	16/32
No fracture	1,648/1,534
Average age at the time of the fracture (years)	4.05 (0-6)
Average BMI at 7 years of agea
Fracture	16.6 (13.1-20.8)
No fracture	16.2 (10.4-34.8)
Location of fracture
Skull	1
Facial bones	1
Vertebral column	1
Clavicle	2
Humerus	10
Radius or ulna	14
One or more phalanges of the hand	4
Other or unspecified parts of femur	4
Tibia or fibula	12

[i] ^aIn the cohort, the body height and weight were recorded when the subjects turned 7 years old. Values are expressed as the average (range) or n. BMI, body mass index.

Genetics of fractures

The GWAS analysis identified one locus with a significant association and six loci with a suggestive result (Table III; Fig. 2). The Q-Q plot of the GWAS is presented in Fig. 3.

Figure 2 Manhattan plot of the genome-wide association study. The red line refers to the level of statistical significance (P=5x10-8) and blue line marks the level of suggestive significance.

Figure 3 Quantile-quantile plot from the genome-wide association study analysis. The red line represents the case where the compared distributions are similar.

Table III Genetic variants identified in the genome-wide association study for childhood fractures performed in the Northern Finland Birth Cohort 1986.

Table III

Genetic variants identified in the genome-wide association study for childhood fractures performed in the Northern Finland Birth Cohort 1986.

SNP ID	Chr	Position	EA	NEA	EAF	BETA	SE	P-value	Adjacent gene(s)
rs112635931	10	11891165	A	G	0.0512	3.187	0.551	7.28x10-9	PROSER2, PROSER2-AS1
rs9827298	3	4023053	G	C	0.0635	2.325	0.437	1.01x10-7	LRRN1, SETMAR, SUMF1
rs374077976	8	126506632	C	CTT	0.1462	1.752	0.333	1.46x10-7	TRIB1, NSMCE2
rs41316954	9	101767382	T	C	0.0890	1.809	0.375	1.38x10-6	COL15A1, TGFBR1
rs17762577	14	97517062	T	C	0.4670	-0.955	0.210	5.53x10-6	VRK1
rs35417231	7	131772615	T	C	0.0703	1.988	0.412	1.37x10-6	PLXNA4
rs111299584	17	69935057	A	AT	0.1007	1.715	0.370	3.50x10-6	SOX9, SOX9-AS1

[i] Chr, chromosome; EA, effect allele; NEA, non-effect allele; EAF, effect allele frequency; SE, standard error; SNP, single-nucleotide polymorphism; PROSER2-AS1, proline- and serine-rich 2 antisense RNA 1; LRRN1, leucine-rich repeat neuronal 1; SETMAR, SET domain and mariner transposase fusion gene; SUMF1, sulfatase-modifying factor 1; TRIB1, tribbles pseudokinas 1; NSMCE2, NSE2 (MMS21) homolog, SMC5-SMC6 complex SUMO ligase; COL15A1, collagen type XV alpha 1 chain; TGFBR1, transforming growth factor beta receptor 1; VRK1, vaccinia-related kinase 1; PLXNA4, Plexin A4; SOX9, SRY-Box transcription factor 9; SOX9-AS1, SOX9 antisense RNA 1.

The lead SNP rs112635931 with the most significant GWA was located within the proline- and serine-rich 2 (PROSER2) and PROSER2-antisense RNA 1 (AS1) genes (Fig. 4). The lead SNP (rs9827298) with results suggestive of an association with fracture risk is located in the intergenic region near the leucine-rich repeat neuronal 1 (LRRN1), SET domain and mariner transposase fusion gene (SETMAR) and sulfatase-modifying factor 1 (SUMF1) genes (Fig. 5). According to GWASdb and GWAS catalog analyses, no genome-wide significant associations with any other trait were reported for either of the SNPs. They were also not present in the RegulomeDB. In addition, the present study did not identify any SNPs associated with other phenotypes that may be in linkage disequilibrium with the lead SNP.

Figure 4 LocusZoom of the rs112635931, Chr 10:11891165, P=7.28x10-9. Each dot represents the P-value of the SNP obtained from the genome-wide association study. The color of the dot represents the linkage disequilibrium (r2), which means that the frequency of association of the different alleles is higher than what would be expected if the loci were independent and associated randomly. Chr, chromosome; SNP, single-nucleotide polymorphism; PROSER2-AS1, proline- and serine-rich 2-antisense RNA 1; USP6NL, Ubiquitin Specific Peptidase 6 N-terminal Like; ECHDC3, Enoyl-CoA Hydratase Domain Containing 3; UPF2, up-frameshift 2; DHTKD1, Dehydrogenase E1 And Transketolase Domain Containing 1; SEC61A2, Selenocystein 61 Translocon Subunit Alpha 2; NUDT5, Nudix Hydrolase 5, Nudix Hydrolase 5; CDC123, Cell Division Cycle 123.

Figure 5 LocusZoom of the rs9827298. Chr3:4023053, P=1x10-7. Each dot represents the P-value of the SNP obtained in the genome-wide association study. The color of the dots represents the linkage disequilibrium (r2), which means that the frequency of association of the different alleles is higher than what would be expected if the loci were independent and associated randomly. Chr, chromosome; SNP, single-nucleotide polymorphism; LRRN1, leucine-rich repeat neuronal 1; SETMAR, SET domain and mariner transposase fusion gene; SUMF1, sulfatase-modifying factor 1.

The SNP rs374077976 was present in GWASdb and this SNP has been associated with total cholesterol in a European population (40). Furthermore, the SNP rs111299584 in chr 3 is considered to likely affect the binding of transcription factors to DNA according to the RegulomeDB. Furthermore, SNPs rs374077976, rs17762577, rs41316954 and rs35417231 were categorized to have minimal evidence of DNA binding of transcription factors in the RegulomeDB. In addition, none of the studied loci were replicated in the summary statistics of the two GEFOS studies (35,38).

Discussion

The results of the present study are important, as the understanding of the genetic associations with childhood fractures is limited, despite the fact that these are common injuries. Recent studies have identified certain risk factors for pediatric fractures both at the micro-level (for example, maternal smoking associates with childhood fractures) (41) and the macro-level (for example, weather conditions associates with upper extremity fractures) (42). However, further studies examining childhood fractures are required (43). Furthermore, previous studies focusing on the possible genetic origins of fractures in pediatric patients without diseases affecting the bone are limited; to the best of our knowledge, the present study provides novel GWAS data on fractures in early childhood.

Previous studies have reported that BMC and BMD are highly heritable (7-11). Furthermore, GWAS have revealed ≥11 different chromosomes that affect BMD in adults (15-21). Chesi et al (8) identified five loci that influence pediatric BMD at multiple skeletal sites. Duren et al (7) identified three genomic regions (2p25.2, 3p25.3 and 17q21.2) with linkage to cortical bone thickness in 10-year-old children (7). These regions contain candidate genes, certain of which have been associated with adult bone mass in other studies (17,44-46). Chesi et al (47) performed a GWAS of areal BMD and BMC and indicated that two loci had a statistically significant association with the BMD and BMC of the distal radius. Furthermore, the bone mass phenotype has been revealed to influence the development of osteoporosis and it is suggested that genetics may impact osteoporotic fracture risk (22).

The primary result of the present study was the statistically significant association between the lead SNP rs112635931 and bone fractures in preschool pediatric patients. Furthermore, SNP rs112635931 was located within the PROSER2 gene, and according to expression quantitative trait loci (eQTL) of the GTEx database, it may regulate the expression of the PROSER2-AS1 gene. Furthermore, the PROSER2 protein is abundant throughout different organs in the human body, including the bone marrow and muscle tissues, based on the GTEx and Human Protein atlas databases (48,49). The expression of PROSER2 protein was also reported to be present in two mouse osteoblast cell lines according to data in the Gene Expression Omnibus database (50,51); however, the exact function of PROSER2 remains to be fully elucidated.

In addition to the associative gene loci, the present study identified six loci that may have an association with fracture risk. The most adjacent genes to the observed locus in chr 3 (lead SNP rs9827298) were LRRN1, SUMF1 and SETMAR. Furthermore, based on the GTEx database, it was speculated that rs9827298 may regulate the expression of LRRN1 in the thyroid, lung and blood; however, the database does not include bone tissue or cells. Furthermore, this SNP may affect the SUMF1 gene, which encodes the formylglycine-generating enzyme (FGE). FGE was detected in the endoplasmic reticulum of bone cells, but also in the brain and skin (33). In addition, FGE converts cysteine into C-α-formylglycine and this activates type I sulfatases (52). Previous studies have indicated that there are ≥8 pathologies caused by disruption of sulfatases, including chondrodysplasia punctata type 1 (53,54). The identified SNP is also in the region of the SETMAR gene, which encodes a fusion protein that binds DNA and functions in DNA repair activities (55). Furthermore, the lead SNP (rs374077976) on chr 8 has previously been associated with total cholesterol (40) and high cholesterol levels have also been suggested to correlate with low BMD in adults (56,57). In addition, there are indications that this may have a relevance in pediatric patients and adolescents (58). The SNP rs41316954 on chr 9 is located within the COL15A1 gene, which encodes the α chain of type XV collagen, a member of the fibril-associated collagens with interrupted helices (FACIT) collagen family. Furthermore, type XV collagen has been revealed to be secreted by human osteoblasts (59). SNP rs111299584 on chr 17 is located upstream of the SOX9 gene, a key transcription factor involved in chondrocyte differentiation. In addition, it has been reported that SOX9 deficiency leads to campomelic dysplasia (60).

A total of two out of the six suggestive loci that were indicated in the GWAS are located in the vicinity of genes that do not have a definitive connection to bone. The first, SNP rs17762577 (chr 14), is located near the vaccinia-related kinase (VRK) serine/threonine kinase 1 gene, which encodes the protein VRK1(61). The expression of this gene is increased in actively dividing cells and it is widely expressed in cells throughout the body (61). The second, SNP rs35417231 on chr 7, encodes a protein (Plexin-A4, a receptor protein) that is necessary for signaling of semaphorins (62).

The major strengths of the present study are the large size of the birth cohort and the valid fracture data from the NHDR. Furthermore, the laboratory and statistical methods used to analyze the potential association were validated and of a high standard. Prior to initiation of the present study, it was determined that less severe childhood fractures may not provide reliable information on underlying disorders of bone development. Thus, only fractures treated in hospital were included to provide data with increased reliability regarding the genetic background associated with the risk of fractures. However, the authors acknowledge that in the future, it may be important to perform another study with a population including conventional out-of-hospital treated low-energy fractures and high-energy fractures. Prior to initiation of the study, it was also decided that only fractures among pre-school aged pediatric patients would be included; in later childhood, individual behavioral factors and recreational activities have a greater impact on the individual risk of conventional bone fractures (63,64).

However, there were certain limitations to the present study. For instance, the genetic data were not available for all patients in the cohort, and thus, the study group did not include all the pediatric patients with fractures treated in hospital prior to the age of 7 years. Furthermore, a larger study population with larger amounts of available genetic material would have increased the strength of the GWAS. However, the present sample of 3,230 is still satisfactory for genetic research to assess the preliminary hypothesis of the study. As another limitation, no subgroup analyses were performed, which may have provided further information on differences between sexes or differences between fractures suffered at different ages (63,64). In addition, the number of cases was modest (n=48) and no further sub-group analyses were performed at this stage. Therefore, to validate the present results, they require to be verified in other populations in future studies.

Bone density measurements were not available in the present study. Therefore, it was not possible to perform any analysis regarding BMC or BMD. The fractures that had been caused by a high-energy injury, e.g. a motor vehicle collision, were not excluded, as information on the mechanism of injury was not available for all fractures. The potential effect of the injury energy on each individual fracture, as a confounding factor, remains elusive and the results of the present study should be interpreted with this limitation in mind. Furthermore, the healing outcomes of the fractures were not available and this was another limitation of the study.

An additional limitation of the present study is the lack of information regarding the other potential confounding factors in early childhood, including mental or physical health conditions or abnormal weight gain. Furthermore, cases with diseases clearly affecting the bone, including malignancies of the bone and osteogenesis imperfecta, were removed from the study, but other types, e.g. growth retardations, were not excluded. In future studies, it may be useful to consider mental deficits, including Attention Deficit Hyperactivity Disorder or developmental retardations, which may increase the risk of fracture (65). However, two essential factors affecting the pediatric fracture risk in general, sex and age, were assessed in the present analyses.

The nutritional status of pediatric patients over the duration of the in-hospital stay has not been previously documented. In Finland, in general, if a pediatric patient had their first or second fracture and there is an obvious explanation, including a fall, vitamin D levels or other laboratory parameters are not assessed. In cases where chronic disease is suspected, further investigations are performed as a common policy in Finnish hospitals. Among the patients with fracture in the present study, no such cases were present. It should also be noted that vitamin D supplementation has been provided to Finnish children since the 1940s (66,67). In addition, normal growth and development of all children in Finland are followed up by mandatory and regular nurse or doctor out-of-hospital visits at a maternity clinic from birth until the age of graduation from junior high school (68). Due to these measures, rickets and other illnesses caused by nutritional deficiencies are rare (66,67). In future studies, assessment of the possible implication of the vitamin D levels may be beneficial, as variations in these levels may be a confounding factor.

In future studies it would be beneficial to also take into account recurrent fractures. In the present study only one patient had suffered from >one fracture in the first six years of life and therefore the recurrence of the fractures was not taken into account in the analysis.

To the best of our knowledge, the present study was the first to evaluate the potential association between pediatric fractures and genetics in this particular setting. The present results suggested that genes may have a role in the fracture risk of children, even in cases without a disease affecting the bones. However, further research in this field is required and future studies should be performed with an additional, larger dataset. In addition, further investigations are required to identify the underlying biological mechanisms of the reported association.

Acknowledgements

The authors would like to acknowledge the biostatistician Mrs. Eeva Vaaramo, M.Sc. (University of Oulu, Faculty of Medicine, Northern Finland Birth Cohort Studies), who was involved in planning the study, and gathering and analyzing the data. The authors would also like to thank Dr Fiona Hanlon-Dearman (University of Alberta, Edmonton AB, Canada) who performed the proofreading of the article.

Funding

No funding was received.

Availability of data and materials

The data that support the findings of this study are available from the NFBC (www.oulu.fi/nfbc/), but restrictions apply to the availability of these data, which were used under license for the present study, and so are not publicly available. The remaining datasets used and/or analzyed are available from the corresponding author on reasonable request and with permission of the NFBC.

Authors' contributions

RP conceived the and designed the present study, and drafted and revised the manuscript. SS acquired, analyses and interpretated the data. SS was also involved in drafting and revising the manuscript. LK was involved in the conception and design of the study, and also took part in the drafting and revising the manuscript. WS was involved in the conception and design of the study as well as revising the manuscript critically. MM was involved in the conception and design of the study. MM also revised the manuscript critically. JJS responsible for the conception and design of the study, and JJS was also involved in drafting and revising the manuscript. All authors read and approved the final manuscript

Ethics approval and consent to participate

The Ethics Committee of the Northern Ostrobothnia Hospital District approved the study (approval no. 12/2003; 27th of Feb 2003). Participation was voluntary and all participants signed a written informed consent form. Personal information was replaced with an identification code and the data were handled at a group level only. The study followed the principles of the Declaration of Helsinki (Oulu, Finland).

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References