Gut microbiota is associated with bone mineral density : an observational and genome-wide environmental interaction analysis in the UK Biobank cohort.

AIMS: Despite the interest in the association of gut microbiota with bone health, limited population-based studies of gut microbiota and bone mineral density (BMD) have been made. Our aim is to explore the possible association between gut microbiota and BMD.
METHODS: A total of 3,321 independent loci of gut microbiota were used to calculate the individual polygenic risk score (PRS) for 114 gut microbiota-related traits. The individual genotype data were obtained from UK Biobank cohort. Linear regressions were then conducted to evaluate the possible association of gut microbiota with L1-L4 BMD (n = 4,070), total BMD (n = 4,056), and femur total BMD (n = 4,054), respectively. PLINK 2.0 was used to detect the single-nucleotide polymorphism (SNP) × gut microbiota interaction effect on the risks of L1-L4 BMD, total BMD, and femur total BMD, respectively.
RESULTS: We detected five, three, and seven candidate gut microbiota-related traits for L1-L4 BMD, total BMD, and femur BMD, respectively, such as genus Dialister (p = 0.004) for L1-L4 BMD, and genus Eisenbergiella (p = 0.046) for total BMD. We also detected two common gut microbiota-related traits shared by L1-L4 BMD, total BMD, and femur total BMD, including genus Escherichia Shigella and genus Lactococcus. Interaction analysis of BMD detected several genes that interacted with gut microbiota, such as phospholipase D1 (PLD1) and endomucin (EMCN) interacting with genus Dialister in total BMD, and COL12A1 and Discs Large MAGUK Scaffold Protein 2 (DLG2) interacting with genus Lactococcus in femur BMD.
CONCLUSION: Our results suggest associations between gut microbiota and BMD, which will be helpful to further explore the regulation mechanism and intervention gut microbiota of BMD. Cite this article: Bone Joint Res 2021;10(11):734-741.

To investigate the possible association of gut microbiota with bone mineral density (BMD).

Genus Escherichia Shigella and genus Lactococcus were associated with L1-L4 BMD, total BMD, and femur total BMD.

Genus Dialister was significantly associated with total BMD via interaction with phospholipase D1 (PLD1) and endomucin (EMCN).

Genus Lactococcus was significantly associated with femur BMD via interaction with Collagen Type XII Alpha 1 Chain (COL12A1) and Discs Large MAGUK Scaffold Protein 2 (DLG2).

This is the first population-based association analysis between gut microbiota with L1-L4 BMD, total BMD, and femur total BMD.

The populations from the UK Biobank and the Flemish Gut Flora Project might differ in ethnic and genetic background, and this might cause possible biases for analysis.

Introduction

Osteoporosis is the most common metabolic bone disorder, caused by an interaction of numerous disease susceptibility genes and environmental factors. The characteristics of osteoporosis are low bone quality, increase in bone resorption, and increased fracture risk. The prevalence of osteoporosis ranges from 6.6% to 22.1% in the EU, and the risk of fragility fracture occurring in the remaining lifetime from 50 years old is 50% for women and 20% for men. Osteoporosis is defined clinically through the measurement of bone mineral density (BMD), which remains the single best predictor of fracture. The affecting factor of BMD has been identified to associate with environment and heredity. For example, twin and family studies have shown that 50% to 85% of the variance was genetically determined in BMD.

Gut microbiota is the whole of commensal, symbiotic, and pathogenic microorganisms living in human gut. The alteration of disease-associated gut microbiota is often characterized by a decrease in species richness and proliferation of microbiota taxa. Gut microbiota is associated with alterations of bone metabolism, bone mineral absorption, and immune regulation in osteoporosis. Studies have shown that gut microbiota regulation has a potential effect on bone mineral density and bone health. For example, the interaction between gut microbiota and the host contributes to the maturation of the host’s immune system, which has an important role in bone metabolism. Abundant evidence suggests that gut microbiota can interact with nonintestinal cells, such as immune cells and dendritic cells, to produce small molecules including short-chain fatty acids, indole derivatives, and polyamines. Receptors of some molecules are expressed on immune cells and regulate the differentiation of T effector cells. Gut microbiota can also increase calcium absorption and regulate the production of serotonin in the gut, a molecule that interacts with bone cells and has been suggested to regulate bone mass. In addition, BMD and immune systems are closely correlated. Interleukin-6 (IL-6) plays an important role in immune response and bone metabolism, enhances macrophage activation and antigen presentation, and mediates the action of osteoblasts and osteoclasts through complex mechanisms. Until now, gut microbiota has been well studied in immune response, but the specific role of gut microbiota species on BMD remains unclear.

Genome-wide association study (GWAS) has succeeded in revealing single-nucleotide polymorphisms (SNPs) that contribute to the associated traits of BMD and osteoporosis. Nevertheless, GWAS results show that the effect sizes of individual causal loci are relatively small. To solve this dilemma, researchers proposed the polygenic risk score (PRS), a score reflecting the sum of all known risk loci. PRS has contributed to the genetic architecture of skeletal disease traits by its ability to predict disease status. Complex human diseases were considered to involve the interaction between environmental and lifestyle factors, as well as inherited susceptibility. The genome-wide environmental interaction (GWEI) study aims to describe the interactions between genetic and environmental factors and the effects on human diseases. Although an association between gut microbiota and BMD has been reported in observational studies, the precise association loci and gut microbiota remain undetermined.

In this study, the UK Biobank data were used to calculate individual PRSs for 114 gut microbiota-related traits. Linear regressions were used to analyze the correlation between each gut microbiota-related PRS with L1-L4 BMD (4,070 participants), total BMD (4,056 participants), and femur total BMD (4,054 participants), respectively. Using the calculated gut microbiota-related PRSs as covariates, GWEI analyses were performed to explore the effects of gene-gut microbiota interactions on the development of BMD.

Methods

Definition of bone mineral density in UK Biobank

The UK Biobank recruited about 500,000 participants aged between 40 and 69 years and conducted prospective studies on them from 2006 to 2010. Briefly, the data field of BMD has three UK Biobank categories, including L1-L4 BMD (Data field 23204), total BMD (Data field 23236), and femur total BMD (Data field 23291). Continuous values of BMD measurement were output from the dual-energy X-ray absorptiometry (DXA) system (g/cm2). BMD was measured at the lumbar spine (L1-L4), total left femur, and total body by DXA. Study subjects reported age, sex, height, and weight on a touchscreen questionnaire. Participants with invalid data on the outcome measure or relevant covariates were excluded in this study. After removing the participants without the calculated gut microbiota-related PRS, 4,070 participants for L1-L4 BMD, 4,054 participants for femur total BMD, and 4,056 participants for total BMD were included for association analysis (Table I).

Table I.

Descriptive characteristics for bone mineral density participants.

Characteristic	L1-L4 BMD	Total BMD	Femur total BMD
Participants, n	4,070	4,056	4,054
Sex, male (%)	1,965 (48.28)	1,954 (48.18)	1,962 (48.40)
Mean age, yrs (SD)	56.05 (7.45)	56.03 (7.44)	56.06 (7.44)
Mean height, cm (SD)	169.93 (9.51)	169.91 (9.51)	169.94 (9.50)
Mean weight, kg (SD)	77.39 (15.29)	77.34 (15.25)	77.39 (15.28)

BMD, bone mineral density; SD, standard deviation.

Genotyping, imputation, and quality control in UK Biobank

Genotyping, imputation, and quality control (QC) for 487,409 individuals were performed by the UK Biobank group. Briefly, the UK BiLEVE Axiom array and UK Biobank Axiom array, which share over 95% of their marker content, were used for genotyping. IMPUTE4 was used for imputation in chunks of about 50,000 imputed markers with a 250 kb buffer region. Marker-based QC was performed using statistical tests designed primarily to check for consistency of genotype calling across experimental factors. Sample-based QC was performed using the metrics of missing rate, heterozygosity, and a set of 15,766 high quality markers on the X and Y chromosomes. More information about genotyping, imputation, QC, and physical measurements has been described previously.

GWAS summary data of gut microbiota

A total of 114 gut microbiota-related traits used in this study were derived from a recent publicly available large-scale GWAS of human gut microbiota. Genetic associations between the human gut microbiota and host genetic variation were identified using the faecal 16S ribosomal RNA gene sequences and human host genotype data. Briefly, the sequencing was carried out for Flemish Gut Flora Project (FGFP) individuals on the Illumina HiSeq platform. Classifications with low confidence at the genus level (< 0.8) were organized in the arbitrary taxon ‘unclassified group’. The DADA2 pipeline yielded count data for 499 taxa across five levels of the microbiota phylogeny, from phylum to genus. After estimating the proportion of gut microbiota variation explained by genetic variation among individuals, associations between genetic variants and specific microbial traits (MTs) were identified by fitting linear, logistic, multinomial, and multivariate regressions assuming an additive genetic model. The Human Core Exome v.1.0 and the Human Core Exome v.1.1 arrays were used for genotyping. Allele calling was performed using GenomeStudio v.2.0.4 (Illumina, USA). In total, 509,886 variants and 2,293 individuals were remained after QC. FGFP genotype data were phased using SHAPEIT3 and imputed with IMPUTE4. Copy number variants were called with PennCNV v.1.0.4. Unique CNVs were defined by unique base pair start and stop locations. SNP variation was linkage disequilibrium (LD)-pruned using PLINK 2.0 and the flag --indep-pairwise 50 5 0.45. In total, 3,321 LD independent loci associated with 16 S gut microbiota species phenotypes were identified to achieve genome-wide significance in the FGFP. Based on 3,321 gut microbiota-related SNPs, 114 gut microbiota-related traits were obtained after eliminating the ones without corresponding SNPs and removing the repetitive gut microbiota-related traits. The detailed information of sample collection, sequencing, microbiome trait preparation, observational analysis, genotyping, heritability, and association analysis are described elsewhere.

Gut microbiota-related PRS calculation and association analysis

According to the standard approach, PLINK 2.0 was used to calculate gut microbiota-related PRS of each study subject using individual genotype data of UK Biobank. Briefly, we set PRS to denote the PRS value of gut microbiota species for the nth subject, defined as:

where l denotes the total number of gut microbiota-associated SNPs; E denotes the effect size of significant gut microbiota-associated SNP i; and D denotes the dosage of the risk allele of the ith SNP for the nth individual (0 is coded for homozygous protective genotype, 1 for heterozygous, and 2 for homozygous polymorphic genotypes). R software (R Foundation for Statistical Computing, Austria) was used to establish linear regression model to evaluate the possible associations between each gut microbiota PRS and target traits of BMD. The PRSs of gut microbiota were set as instrumental variables, while age, sex, height, and weight were set as covariates.

Statistical analysis

The genotype data of BMD were firstly adjusted for age and sex, and ten principal components of population structure (PCs) using linear regression models, and the residuals from the regression model were then used for GWEI analysis, respectively. The command ‘glm’ of PLINK 2.0 was used to analyze the interaction between SNPs and the PRS of significant gut microbiota for BMD, setting PRSs as covariates. For quality control, we removed the SNPs with call rates < 90%, Hardy-Weinberg equilibrium (HWE) < 0.001, or minor allele frequencies (MAF) < 0.01. The kinship coefficients were estimated by KING software (University of Virginia, USA) to remove the genetically related subjects. Rectangular Manhattan plot and QQ plot were produced using the “CMplot” package (https://github.com/YinLiLin/R-CMplot) in R platform.

Results

Associations of gut microbiota with bone mineral density

Linear regression detected five candidate gut microbiota associated with L1-L4 BMD in UK Biobank, such as genus Dialister_HB (p = 0.004), genus Dialister_RNT (p = 0.005), and genus Escherichia Shigella_HB (p = 0.007) (Table II). For total BMD, we detected three candidate gut microbiota, including genus Lactococcus_HB (p = 0.013), genus Escherichia Shigella_HB (p = 0.035), and genus Eisenbergiella_HB (p = 0.046) (Table III). For femur total BMD, we detected seven candidate gut microbiota, such as genus Escherichia Shigella_HB (p = 0.005), genus Lactococcus_HB (p = 0.021), and genus Dialister_RNT (p = 0.028) (Table IV).

Table II.

The gut microbiota associated with L1-L4 bone mineral density.

Gut microbiota	β	SE	p-value*
G Dialister_HB	0.203	0.070	0.004
G Dialister_RNT	-0.170	0.061	0.005
G Escherichia Shigella_HB	-0.073	0.027	0.007
G Lactococcus_HB	-0.195	0.076	0.010
G Senegalimassilia_HB	-0.072	0.035	0.039

Linear regression.

G, genus; HB, hurdle binary; RNT, rank-normal transformation; SE, standard error.

Table III.

The gut microbiota associated with total bone mineral density.

Gut microbiota	β	SE	p-value*
G Lactococcus_HB	-0.117	0.047	0.013
G Escherichia Shigella_HB	-0.035	0.017	0.035
G Eisenbergiella_HB	-0.110	0.055	0.046

Linear regression.

G, genus; HB, hurdle binary; SE, standard error.

Table IV.

The gut microbiota associated with femur total bone mineral density.

Gut microbiota	β	SE	p-value*
G Escherichia Shigella_HB	-0.059	0.021	0.005
G Lactococcus_HB	-0.136	0.059	0.021
G Dialister_RNT	-0.104	0.047	0.028
G Dialister_HB	0.119	0.054	0.029
G Veillonella_HB	-0.107	0.049	0.029
F Enterobacteriaceae_HB	0.064	0.031	0.037
F Veillonellaceae_HB	-0.091	0.044	0.039

Linear regression.

F, family; G, genus; HB, hurdle binary; RNT, rank-normal transformation; SE, standard error.

We further compared the above association analysis results, and found two candidate gut microbiota shared by L1-L4 BMD, total BMD, and femur total BMD, including genus Escherichia Shigella_HB (P L1-L4 BMD = 0.007, P total BMD = 0.035, P femur total BMD = 0.005), and genus Lactococcus_HB (P L1-L4 BMD = 0.010, P total BMD = 0.013, P femur total BMD = 0.021). In addition, two candidate gut microbiota were shared by L1-L4 BMD and femur total BMD, including genus Dialister_HB (P L1-L4 BMD = 0.004, P femur total BMD = 0.029), and genus Dialister_RNT (P L1-L4 BMD = 0.005, P femur total BMD = 0.028).

Interaction analysis of gut microbiota with bone mineral density

For L1-L4 BMD, we detected one significant SNP interacted with genus Escherichia Shigella_HB, rs74862545 (CCND2-AS1, p = 1.65 × 10-8) (Figure 1a and 1b). For total BMD, we detected 19 significant SNPs interacted with genus Dialister_RNT, such as rs79540008 (PLD1, p = 3.09 × 10-8), rs78658424 (endomucin (EMCN), p = 1.62 × 10-9), and rs191380733 (HSPA7, p = 1.07 × 10-8) (Figure 1c and 1d). The significant GWEI results (p < 5.00 × 10−8) of total BMD are summarized in Supplementary Table i. For femur total BMD, we detected 76 significant SNPs interacted with genus Lactococcus_HB, such as rs191860862 (COL12A1, p = 1.44 × 10-11), rs74777764 (Discs Large MAGUK Scaffold Protein 2 (DLG2), p = 2.90 × 10-8), and rs111824870 (nuclear factor of activated T cells C2 (NFATC2), p = 4.84 × 10-8) (Figure 1e and 1f). The significant GWEI results (p < 5.00 × 10−8) of femur total BMD are summarized in Supplementary Table ii. We also detected one significant SNP interacted with genus Escherichia Shigella_HB, rs190533440 (p = 1.84 × 10-8) (Figure 1g and 1h).

Fig. 1

Genome-wide environmental interaction (GWEI) study of bone mineral density. a) and b): GWEI study in genus Escherichia Shigella_HB of L1-L4 bone mineral density (BMD). c) and d): GWEI study in genus Dialister_RNT of total BMD. e) and f): GWEI study in genus Lactococcus_HB of femur total BMD. g) and h): GWEI study in genus Escherichia_Shigella_HB of femur total BMD. In Manhattan plot, the black solid line indicates the p-value threshold for genome-wide significance (p < 5 × 10−8), while the black dotted line indicates p-value threshold for suggestive significance (p < 5 × 10−5). In QQ plot, a graphical representation of the deviation of the observed p-values from the null hypothesis: the observed p-values for each single-nucleotide polymorphism (SNP) are sorted from largest to smallest and plotted against expected values from a theoretical χ2-distribution. HB, hurdle binary; RNT, rank-normal transformation.

Discussion

In this study, a recent large-scale GWAS was used to obtain gut microbiota-associated loci. The UK Biobank data were used to conduct PRS analysis of L1-L4 BMD, total BMD, and femur total BMD for each individual in the UK Biobank cohort, respectively. The GWEI analyses were performed to detect candidate SNP × gut microbiota interaction effects on L1-L4 BMD, total BMD, and femur total BMD, respectively. Our study observed associations of gut microbiota with BMD, and detected several candidate loci that interacted with gut microbiota for BMD.

Measurements of DXA bone mass are considered to be the “gold standard” for the diagnosis of osteoporosis, but subtle differences exist when different bone sites were measured. DXA can be used to measure BMD at multiple sites. The most frequently used sites include L1-L4 BMD, total BMD, and/or femur total BMD. Several studies have explored the best location for DXA scans to measure BMD and diagnose osteoporosis. Compared with total DXA, the lumbar spine DXA was associated with a higher prevalence of osteoporosis. Additionally, the International Society for Clinical Densitometry advocates the measurement of BMD in the first four lumbar spines, and the diagnosis of osteoporosis based on the lowest T-score at the measurement locations. However, since BMD varies by age and site, no single measurement point can detect BMD in all cases. For example, the rate of bone loss is influenced by the patient’s age and bone sites. During the perimenopausal and early post-menopausal period, bone loss occurs primarily in the spine due to the effects of oestrogen deficiency on the trabecular bone. Thus, measurement of the hip alone may miss the diagnosis of osteoporosis in this group of patients. Conversely, in older adults, structural changes in the back of the spine, such as vascular calcification and degenerative arthritis, may mistakenly increase bone density in the spine and thereby limit its utility. To date, no study has been performed on the relationship between BMD and gut microbiota in different bone sites. In this study, we selected three common BMD phenotypes to explore the relationship between them and gut microbiota.

The genus Escherichia Shigella and genus Lactococcus were observed to associate with L1-L4 BMD, total BMD, and femur total BMD in this study. Recently, it has been demonstrated that the interactions between gut microbiota and phosphorus/vitamin D were crucial in bone development and mineralization. The positive effects of Escherichia Faecium on phosphorus metabolism were associated with changes in Escherichia Shigella in broilers, which improved gut phosphorus absorption and bone-forming metabolic activities, and decreased phosphorus excretion. However, bone minerals and apatite also serve as a dumping ground for trace elements and drugs, which seriously affects the bone health. Vitamin D3 supplementation decreased the relative abundance of Escherichia Shigella in the upper gastrointestinal tract. Further, Yang et al found that Shigella flexneri could induce robust inflammasome activation in mouse bone marrow macrophages. The genus Escherichia Shigella may have a vital role in BMD by affecting bone development and mineralization.

The genus Lactococcus comprises 12 species. Kimoto-Nira et al found that genus Lactococcus was related to BMD changes. Oral administration of heat-killed Lactococcus Lactis to aged SAMP6 mice (a senescence-accelerated mice strain that develops osteoporosis with ageing) was associated with reduced bone density loss. Similarly, the abundance of genus Lactococcus was markedly decreased in mice with osteoporosis. Shimada et al further reported the isolation of an enzyme related to daidzein metabolism and equol production in Lactococcus strain, which is more potent than that of other isoflavones on ameliorative ability against lower BMD. Considering the growing evidence demonstrating that gut microbiota is related to osteoporosis, our results indicate that genus Escherichia Shigella and genus Lactococcus may associate with the alteration of BMD.

Interaction analysis of total BMD detected that phospholipase D1 (PLD1) had interaction effects with genus Dialister. Phospholipases are suspected to affect bone remodelling and formation, as evidenced by their expression and activity in forming osteoblasts and chondrocytes, and resorbing osteoclasts. The different isoforms of PLD in chondrocytes and osteoblasts were previously reported to regulate differentiation, maturation, and function of cells. For instance, Yoo et al suggested that the targeting inhibition of PLD1 could ameliorate bone erosion and cartilage destruction by suppressing osteoclastogenesis. EMCN is another candidate gene that had interaction effects with genus Dialister in total BMD. Osteogenesis during bone modelling and remodelling is coupled with angiogenesis. EMCN interferes with the assembly of focal adhesion complexes and inhibits interaction between cells and the extracellular matrix. Although there is less evidence to link genus Dialister to the regulation of EMCN and PLD1, our results suggest that genus Dialister may influence the total BMD by affecting the expression of EMCN and PLD1. Our interaction analysis of femur total BMD highlighted COL12A1 as one of the significant genes interacting with genus Lactococcus. Collagen XII is the largest member of fibril-associated collagens with interrupted triple helix family, assembled from three identical α-chains encoded by the COL12A1 gene. Bone formation is precisely regulated by cell-to-cell communication in osteoblasts. COL12A1 was found to downregulate in aged osteoblasts. Izu et al demonstrated that genetic deletion of COL12A1 impaired osteoblast connection and/or communication in mice, resulted in reduced bone mass, and increased bone fragility.

DLG2 and NFATC2 were also identified to have interaction effects with genus Lactococcus in femur total BMD. DLG2 was found to be associated with the new bone formation in RNA sequencing. NFATC2 is important for the immune response, whereas NFATC1 is a crucial transcription factor for osteoclast differentiation and osteoclastogenesis in vitro. Bone formation was inhibited in NFATC1- and NFATC2-deficient cells, and stimulated in NFATC1 overexpression cells, suggesting that NFATC1 and NFATC2 were associated with osteoblastic bone formation and osteoporosis. Additionally, NFATC2 activation in osteoblasts could inhibit bone formation and cause cancellous bone osteopenia. While the considerable associations of genus Lactococcus with SNPs were reported in BMD, the causal relationships and biological mechanisms remain elusive. From a genetic perspective, our research suggests a possible effect of genus Lactococcus on femur total BMD.

There are several limitations in this study. Firstly, given the lack of information about culture/geographical background, a measure of income, education, or socioeconomic status, as well as dietary habit, we could not consider these potential confounding factors in our analysis. Secondly, gut microbiota data were taken from the FGFP, while BMD data were taken from the UK Biobank. The tiny differences in demographic backgrounds may partly skew our results. Thirdly, although the gut microbiota and GWEI reported in this study are significantly related to L1-L4 BMD, total BMD, and femur total BMD, which is consistent with some previous evidence, further experimental studies are needed to explore and confirm the underlying molecular biological mechanisms. Finally, the GWAS and gut microbiota data in this study were obtained from individuals of European ancestry, which should be applied to other races with care.

In summary, we performed PRS and GWEI analysis to evaluate the associations between gut microbiota and L1-L4 BMD, total BMD, and femur total BMD. Our findings suggest the potential role of gut microbiota in the aetiology of osteoporosis and, in particular, we identified that genus Dialister and genus Lactococcus may have significant effects on total BMD and femur BMD, respectively. This also highlights that many loci and genes involved in gut microbiota-BMD interactions are yet to be characterized, and that studying gut microbiota in the context of their interactions with other diseases offers a way to discover new areas of medicine. Furthermore, it is suggested by this study that gut microbiota might be a therapeutical target of bone diseases.